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Abstract 

Three-dimensional models which contain both geometry and texture have numerous appli- 
cations such as urban planning, physical simulation, and virtual environments. A major fo- 
cus of computer vision (and recently graphics) research is the automatic recovery of three- 
dimensional models from two-dimensional images. After many years of research this goal is 
yet to be achieved. Most practical modeling systems require substantial human input and 
unlike automatic systems are not scalable. 

This thesis presents a novel method for automatically recovering dense surface patches us- 
ing large sets (1000's) of calibrated images taken from arbitrary positions within the scene. 
Physical instruments, such as Global Positioning System (GPS), inertial sensors, and incli- 
nometers, are used to estimate the position and orientation of each image. Essentially, the 
problem is to find corresponding points in each of the images. Once a correspondence has been 
established, calculating its three-dimensional position is simply a matter of geometry. Long 
baseline images improve the accuracy. Short baseline images and the large number of images 
greatly simplifies the correspondence problem. The initial stage of the algorithm is completely 
local and scales linearly with the number of images. Subsequent stages are global in nature, 
exploit geometric constraints, and scale quadratically with the complexity of the underlying 
scene. 

We describe techniques for: 1) detecting and localizing surface patches; 2) refining camera 
calibration estimates and rejecting false positive surf els; and 3) grouping surface patches into 
surfaces and growing the surface along a two-dimensional manifold. We also discuss a method 
for producing high quality, textured three-dimensional models from these surfaces. Someofthe 
most important characteristics of this approach are that it: 1) uses and refines noisy calibra- 
tion estimates; 2) compensates for large variations in illumination; 3) tolerates significant soft 
occlusion (eg. tree branches); and 4) associates, at a fundamental level, an estimated normal 
(eliminating the frontal -planar assumption) and texture with each surface patch. 
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Chapter 1 
Introduction 



Three-dimensional models which contain both geometry and texture have nu- 
merous applications. They a re used extensively as virtu a I environments for en- 
tertainment (e.g. games, videos, commercials and movies). Three-dimensional 
models also have serious applications such as physical simulation and urban 
planning. Recent advances in hardware and software have improved our abil- 
ity to store, navigate and display large, complex models. However, model con- 
struction is frequently a tedious and time consuming task requiring significant 
human input. As the size, complexity, and realism of models increase, man- 
ual, and even interactive or semi-automatic, model construction will rapidly 
become impractical. 

This thesis presents a novel method for automatically recovering dense sur- 
face patches usi ng I arge sets (1000's) of cal i brated i mages taken from arbitrary 
positions within the scene. Physical instruments, such as Global Positioning 
System (GPS), inertial sensors, and inclinometers, are used to estimate the po- 
sition and orientation of each image. The large set of images acquired from a 
wi de range of vi ewpoi nts ai ds i n i dentifyi ng correspondences and enabl es accu- 
rate computation of their three-dimensional position. We describe techniques 
for: 

• Detecting and localizing surface patches. 

• Refining camera calibration estimates and rejecting false positive surfels. 

• Grouping surface patches into surfaces. 

• Growing surfaces along a two-dimensional manifold. 

In addition, we also discuss a method for producing high quality, textured 
three-dimensional models from these surfaces. The initial stage of the al- 
gorithm is completely local making it easily parallelizable. Some of our ap- 
proach's most important characteristics are: 

• It is fully automatic. 

• It uses and refines noisy calibration estimates. 
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CHAPTER 1. INTRODUCTION 



It compensates for large variations in illumination. 

It matches image data directly in three-dimensional space. 

It tolerates significant soft occlusion (eg. tree branches). 

It associates, at a fundamental level, an estimated normal (eliminating 
thefrontal-planar assumption) and texture with each surface patch. 




Figure 1-1: Albrecht Du'rer, Artist Drawing a Lute 1525. 



1.1 Background 



The basics of perspective image formation have been known for more than 2000 
years and date back to Pythagoras, Euclid, and Ptolemy. In the 15 th century, 
LeonoAlberti published thefirst treatise on perspective and later that century 
Leonardo da Vinci studied perspective projection and depth perception. By the 
16 th century these concepts were well known to artists (eg. Figurel-1). Inthe 
18 th centuryj ohan Heinrich Lambert developed a technique called spaceresec- 
tion to find the point in space from which a picture was made. Through the end 
of the 18 th century, the study of perspective focused exclusively on image for- 
mation. The main motivation was accurately depicting the three-dimensional 
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world on a two-dimensional canvas. Early in the 19 th century the photograph 
was invented and the natural question followed: Can two-dimensional pho- 
tographs be used to deduce thethree-dimensional world which produced them? 
The problem of recovering three-dimensional information from a set of pho- 
tographs or images is essentially the correspondence problem: Given a point in 
one image find the corresponding point in each of the other images. Typically, 
photogrammetric approaches (Section 1.1.1) require manual identification of 
correspondences, while computer vision approaches (Section 1.1.2) rely on au- 
tomatic identification of correspondences. If the images are from nearby po- 
sitions and similar orientations (short baseline), they often vary only slightly, 
simplifying the identification of correspondences. Once sufficient correspon- 
dences have been identified, solving for the depth is simply a matter of geom- 
etry. This applies to both calibrated and uncalibrated images. For calibrated 
images (known internal calibration, camera position, and orientation) a sin- 
gle correspondence is sufficient to triangulate the three-dimensional position 
of the scene point which gave rise to the corresponding image points. The un- 
calibrated case is more complicated requiring additional correspondences to 
recover the calibration as well as the three-dimensional information. 




Figure 1-2: Calculating Three-Dimensional Position. 

Real images are composed of noisy, discrete samples, therefore the calcu- 
lated three-dimensional location of a correspondence will contain error. This 
error is a function of the baseline or distance between the images. Figure 1- 
2 shows how the location of point P can be calculated given two images taken 
from known cameras Ci and C 2 and corresponding points p 1 andp 2 within those 
images, which are projections of P. The exact location of p 1 in the image is un- 
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certain, as a result P can lie anywhere within the left cone. A similar situation 
exists for p 2 . If p 1 and p 2 are corresponding points, then P could lie anywhere 
in the shaded region. Clearly, for this situation increasing the baseline be- 
tween Ci and C 2 will reduce the uncertainty in the location of P. This leads 
to a conflict: short baselines simplify the identification of correspondences but 
produce inaccurate results; long baselines produce accurate results but compli- 
cate the identification of correspondences. Real images also frequently contain 
occlusion and significant variations in illumination both of which complicate 
the matching process. 

1.1.1 Photogrammetry 




Figure 1-3: 19 th Century Stereoscope. 

About 1840, Dominique Francois J ean Arago, the French geodesist, advo- 
cated the use of photography by topographers. A decade I ater, Ai me Laussedat, 
a Colonel in the French Army Corps of Engineers, set out to prove that photog- 
raphy could be used to prepare topographic maps. Over the next 50 years, his 
work was so complete that Aime Laussedat is often referred to as the father 
of photogrammetry. Prior to the advent of computing, analog techniques were 
used to physically reprqject single images and stereo pairs. The stereoscopic 
viewer, shown in Figure 1-3, is one such reprqjection device. The availability 
of computing in the mid 20 th century enabled the introduction of analytical 
techniques. For the first time multiple (more than two) photographs could be 
analyzed simultaneously [Wolf, 1974, Slama etal., 1980]. 

In spite of the fact that it originated with ground based imagery, modern 
photogrammetry uses long-range aerial imagery almost exclusively. This is 
in contrast with the close-range ground base imagery used as input for this 
thesis. For accurate results, photogrammetry requires good camera calibra- 
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tion. The internal parameters 1 are measured in a controlled environment 
and traditionally the external parameters are estimated using points whose 
three-dimensional position is known or ti e-poi nts. Classical photogrammetry 
requires a large amount of human input, typically in the form of identifying 
corresponding features in multiple photographs. The correspondences are then 
used to recalculate the external parameters, as well as determine the three- 
dimensional position of selected features. Two classes of techniques are used 
to reestimate the external parameters. Closed-form analytic techniques re- 
quire fewer correspondences and do not need an initial estimate, but tend to 
be numerically unstable and sensitive to error in the data. Global optimization 
techniques, such as bundle-adjustment, require both an initial estimate and 
many correspondences, but generally produce more stable results. The cali- 
bration updates described in Section 5.1 fall into the latter category. Recently, 
automated photogrammetric methods have been explored [G reeve, 1996]. 

A number of the recent interactive modeling systems are based upon pho- 
togrammetry 2 . Research projects such as RADIUS [Collins et al., 1995] and 
commercial systems such as FotoG [Vexed, 1997] are commonly used to ex- 
tract three-dimensional models from images. Good results have been achieved 
with these systems, however the requirement for human input limits the size 
and complexity of the recovered model. One approach to reducing the amount 
of human input is to exploit geometric constraints. The geometric structure 
typical of urban environments can be used to constrain the modeling process. 
Becker and Bove[1995] manually identify groups of parallel and perpendicular 
lines across small sets of images. Shum et al.[1998] interactively draw points, 
lines and planes onto a few panoramic mosaics. Debevecet al.'s Facade system 
[1996] uses a small set of building blocks (cube, cylinder, etc.) to limit the set 
of possible models from a small set of images (at most a couple dozen). The 
major strength of these systems is the textured three-dimensional model pro- 
duced. Debevec et al. cite a hundred fold decrease in human input using the 
block-based approach. Despite this reduction, each image must be processed 
individually by a human to produce a three-dimensional model, making it dif- 
ficult to extend these systems to large sets of i mages. 

1.1.2 Computer Vision 

The field of computer vision began with the advent of computing in the mid- 
20 th century and in many ways developed in parallel to, but separate from 
photogrammetry [Horn, 1986, Mayhew and Frisby, 1991, Faugeras, 1993]. A 
major focus of computer vision is the automatic recovery of three-dimensional 



1 See Appendix A for a discussion of internal and external parameters. 
2 Some might consider these to be in the field of computer vision. We have placed them in 
this section because, liketraditional photogrammetric methods, they require human input. 



18 CHAPTER 1. INTRODUCTION 

information from two-dimensional images. Both calibrated and uncalibrated 
images can be used and as noted above, recovering three-dimensional informa- 
tion can be reduced to finding correspondences 

U ncal i brated I magery 

Several researchers [Longuet-Higgins, 1981, Horn, 1987, Mohr and Arbogast, 
1991, Hartley, 1992, Hartley et al., 1992] have shown that the relative cam- 
era orientations 3 and the three-dimensional structure of the scene can be re- 
covered from images with known internal calibration and unknown external 
calibration. If the internal calibration is also unknown, Faugeras [1992] has 
shown that the scene structure can only be determined up to an unknown pro- 
jective transformation. Structure from motion and direct motion estimation 
are two popular approaches which use uncalibrated imagery. Structure from 
motion techniques [Ullman, 1978, Ullman, 1979, Tomasi and Kanade, 1992, 
Azarbayejani and Pentland, 1995, Seales and Faugeras, 1995] identify and 
track high-level features (eg. edges) across several images. Because both 
structure and motion are recovered, a large set of features must be tracked 
for a robust solution. Direct motion estimation techniques [Mohr et al., 1993, 
Szeliski and Kang, 1994, Ayer and Sawhney, 1995, Adelson and Weiss, 1996, 
I rani et al., 1997] use the optical flow of each pixel across images to directly 
estimate structure and motion. Computing optical flow tends to be sensitive 
to changes in illumination and occlusion. The need to track features or calcu- 
late optical flow across images implies that adjacent images must be close. The 
cost of acquiring such an image set rapidly becomes prohibitive for large-scale 
models. 

C al i brated I magery 

A popular set of approaches for automatically finding correspondences from 
calibrated images is relaxation techniques [Marr and Poggio, 1979]. These 
methods are generally used on a pair of images; start with an educated guess 
for the correspondences; then update them by propagating constraints. These 
techniques often exploit global constraints such as smoothness [Pol lard et al., 
1985] and ordering [Ohta and Kanade, 1985]. They don't always converge and 
don't always recover the correct correspondences. Utilizing only one pair of 
images at a time has several disadvantages. 

• The trade off between easily identifying correspondences and accurate 
results (discussed above) means that these methods generally produce 
less accurate results. 



3 External calibration in an unknown coordinate frame. 
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• The results from one pair are restricted to the regions visible in both im- 
ages. To reconstruct a large volume many pairs are required, leading to 
the difficult problem of integrating the results from multiple pairs. 

Another approach is to use multiple images. Several researchers, such as 
Yachida [1986], have proposed trinocular stereo algorithms and some have 
proposed using a third image to check correspondences [Faugeras and Robert, 
1994]. Others have also used special camera configurations to aid in the cor- 
respondence problem, [Tsai, 1983, Bolles et al., 1987, Okutomi and Kanade, 
1993]. Bolles etal. [1987] proposed constructing an epi polar-plane image from 
a large number of images. In some cases, analyzing the epi polar-plane image 
is much simpler than analyzing the original set of images. The epi polar-plane 
image, however, isonly defined for a limited set of camera positions. Baker and 
Bolles [1989] have also attempted to extend the limited set of usable camera 
positions. Tsai [1983] and Okutomi and Kanade [Okutomi and Kanade, 1993, 
Kanade and Okutomi, 1994] defined a cost function which was applied directly 
to a set of images. The extremum of this cost function was then taken as the 
correct correspondence. Occlusion is assumed to be negligible. I n fact, Okutomi 
and Kanade state that they "invariably obtained better results by using rela- 
tively short baselines." This is likely the result of using an i mage space match- 
ing metric (a correlation window) and ignoring perspective distortion and oc- 
clusion. Both methods use small sets of images, typically about ten. They also 
limit camera positions to special configurations. Tsai uses a localized planar 
configuration with parallel optical axes; and Okutomi and Kanade use short 
linear configurations. 

Some recent techniques use multiple images which are not restricted to 
rigid camera configurations. Kang et al. [1995] use structured light to aid in 
identifying correspondences. Kang and Szeliski [1996] track a large number of 
features across a small number of closely spaced (~ 3 inches) panoramic images 
(360° field of view). Neither of these approaches is well suited to extended urban 
environments. Structured lighting is difficult to use outdoors and tracking a 
dense set of features across multiple large datasets is difficult at best. 

Collins [1996] proposed a space-sweep method which performs matching in 
three-dimensional space instead of image space. There are several advantages 
to matching in three-dimensional space. It naturally operates on multiple im- 
ages and can properly handle perspective distortion and occlusion, improving 
the matching process. Col I ins' formulation assumes that there is a singleplane 
which separates all of the cameras from the scene. As the plane is swept 
away from the cameras, features from all of the images are projected onto it. 
Sweep plane locations which have more than a minimum number of features 
are retained. The resulting output is a sparse set of three-dimensional points. 
Collins demonstrates his method on a set of seven aerial photographs. Seitz 
and Dyer [1997] proposed a voxel coloring method which also performs match- 
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ing in three-dimensional space. Voxels are reconstructed in a special order such 
that if point P occludes point Q for any image then P is reconstructed before Q. 
To meet this ordering requirement, no scene points are all owed within thecon- 
vexhull of the cameras. Collins' sweep-plane is a special case of this ordering 
requirement. Neither of these approaches issuitablefor reconstructions using 
images acquired from within the scene. 

Kutulakos and Seitz [1998b, 1998a] have proposed an extension of voxel 
coloring called space carving. The separability requirement is removed. An 
initial estimate of the reconstructed volume is required and must be a rela- 
tively tight superset of theactual volume. Voxels are tested in order similar to 
Seitz and Dyer. Background pixels are removed and matching is performed on 
raw pixel values. Ifthe projections of a voxel into the images are not consistent, 
the voxel is carved away. The reconstruction is complete when no more voxels 
can be removed. The carving operation is brittle. If a single image disagrees 
the voxel is removed, making this approach particularly sensitive to camera 
calibration errors, illumination changes, complex reflectance functions, image 
noise, and temporal variations. I n short, it is not well suited to outdoor urban 
environments. 

While computer vision techniques eliminate the need for human input in 
general they have several characteristics which limit their applicability in cre- 
ating large, complex, and realistic three-dimensional models. 

• They are frequently fragile with respect to occlusion and variations in 
illumination. 

• They frequently operate on only a few images and do not scale to large 
sets of i mages. 

• The output (eg. depth maps and isolated edges) are generally not directly 
useful as a model. 

Recently, the last item has been partially addressed by Fua [Fua, 1995, Fua 
and Leclerc, 1995, Fua and Leclerc, 1996] who starts with a set of three- 
dimensional points and fits small surface elements to the data. Fitting is ac- 
complished by optimizing an image-based objective function. This is in contrast 
with the approach presented in this thesis which recovers surface elements di- 
rectly from the image data. In addition, as presented Fua's approach does not 
have the ability to fill in holes in the three-dimensional points used as input. 

Finally, several researchers have used probability theory to aid in recov- 
ering three-dimensional information. Cox et al. [Cox, 1994, Cox et al., 1996] 
proposed a maximum-likelihood framework for stereo pairs, which they have 
extended to multiple images. This work attempts to explicitly model occlusions, 
although in a somewhat ad hoc manner. Belhumeur [Belhumeur and Mumford, 
1992, Belhumeur, 1993, Belhumeur, 1996] develops a detailed Bayesian model 



1.1. BACKGROUND 21 

of image formation and the structure of the world. These two approaches serve 
as the inspiration for the probabilistic elements of this thesis. 



1.1.3 Discussion 

Photogrammetric methods produce high quality models, but require human 
input. Computer vision techniques function automatically, but generally do 
not produce usable models, operate on small sets of images and frequently 
are fragile with respect to occlusion and changes in illumination. The work 
presented in this thesis draws from both photogrammetry and computer vi- 
sion. Like photogrammetric methods we produce high quality textured models 
and like computer vision techniques our method is fully automatic. The major 
inspiration derives from Bolles et al. [Bolles et al., 1987, Baker and Bolles, 
1989]. We define a construct called an epi polar image and use it to analyze 
evidence about three-dimensional position and orientation. LikeTsai [1983] 
and Okutomi and Kanade [Okutomi and Kanade, 1993, Kanadeand Okutomi, 
1994] we define a cost function that is applied across multiple images, how- 
ever, we do not evaluate the cost function in image space. Instead, I ike Col I ins 
[1996], Seitz and Dyer [1997], and Kutulakos and Seitz [Kutulakos and Seitz, 
1998b, Kutulakos and Seitz, 1998a] we perform matching in three-dimensional 
space. We also model occlusion and the imaging process similar to Cox et 
al. [1996] and Belhumeur [Belhumeur and Mumford, 1992, Belhumeur, 1993, 
Belhumeur, 1996]. 

There are al so several i mportant differences from previ ous work. The epi po- 
lar image we define is valid for arbitrary camera positions within the scene and 
is capable of analyzing very large sets of images. Our focus is recovering built 
geometry (architectural facades) in an urban environment. However, the al- 
gorithms presented are generally applicable to objects that can be modeled by 
small planar patches. Surface patches (geometry and texture) or surfds are 
recovered directly from the image data. In most cases, three-dimensional posi- 
tion and orientation can be recovered using purely local information, avoiding 
the computational costs of global constraints. Some of the significant charac- 
teristics of this approach are: 

• Large sets of images contain both long and short baseline images and ex- 
hibit the benefits of both (accuracy and ease of matching). It also makes 
our method robust to sensor noise and occlusion, and provides the infor- 
mation content required to construct complex models. 

• Each image is calibrated - its position and orientation in a single global 
coordinate system is estimated. The use of a global coordinate system 
allows data to be easily merged and facilitates geometric constraints. 
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• The fundamental unit is a textured surface patch and matching is done in 
three-dimensional space. This eliminates the need for the frontal-planar 
assumption made by many computer vision techniques and provides ro- 
bustness to soft occlusion (eg. tree branches). Also, the surface patches 
are immediately useful as a coarse model and readily lend themselves to 
aggregation to form more refined models. 

• The algorithm tolerates significant noise in the calibration estimates and 
produces updates to those esti mates. 

• The algorithm corrects for changes in illumination. This allows raw image 
properties (eg. pixel values) to be used avoiding the errors and ambigui- 
ties associated with higher level constructs such as edges or corners. 

• The algorithm scales well. The initial stage is completely local and scales 
linearly with the number of images. Subsequent stages are global in na- 
ture, exploit geometric constraints, and scale quadratically with the com- 
plexity of the underlying scene. 4 

As noted above, not all of these characteristics are unique, but their combi- 
nation produces a novel method of automatically recovering three-dimensional 
geometry and texture from large sets of images. 



1.2 City Scanning Project 

This thesis is part of the City Scanning project whose primary focus is the 
Automatic Population of Geospatial Databases (APGD). As its name implies, 
the main target of the City Scanning project is urban environments such as 
shown in Figure 1-4. In addition to the one presented in this thesis, other 
approaches to three-dimensional reconstruction have been explored as part of 
the City Scanning project. Coorg[1998] hypothesizes large (building size) ver- 
tical faces from sparse edge information and Chou [1998] attempts to deduce 
faces by aggregating sparse three-dimensional features. I n contrast with these 
approaches, this thesis performs dense reconstruction directly from the image 
data. Another major focus of the project is the interactive navigation and ren- 
dering of city sized databases. The following sections describe the acquisition 
and preprocessing phases which produce the calibrated images used as input 
for the algorithms presented in this thesis. For a more complete description of 
the project seetTeller, 1998a, Teller, 1998b]. 



4 This is the worst case complexity. With spatial hashing the expected complexity is linear 
in the number of reconstructed surfels. 
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Figure 1-4: Goal of the city project. Rendering from an idealized model built 
with human input. 
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Figure 1-5: Argus. 




Figure 1-6: Hemispherical tiling for a node. 
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1.2.1 Node Acquisition 

A self contained platform called Argus (Figure 1-5) continues to be developed 
to acquire calibrated imagestDe Couto, 1998]. At the center of Argus is a 
precision computer controlled pan-tilt head upon which a digital camera is 
mounted. The camera is initially calibrated using Tsai's method [Tsai, 1983, 
Tsai, 1987] and is configured to rotate about its center of projection. At each 
node images are acquired in a hemispherical tiling similar to that shown in 
Figure 1-6. Argus is also equipped with a number of navigational sensors in- 
cluding GPS, inertial sensors, odometry and inclinometers which enable it to 
estimate the position and orientation of each image. Currently node positions 
and orientations are also surveyed by hand which serves as validation. The 
raw image data, exposure parameters, and absolute position and orientation 
estimates are recorded by an on-board computer. 




Figure 1-7: Spherical Mosaic for a node. 



1.2.2 Mosaicing and Registration 

The position and orientation estimates obtained during the acquisition phase 
are good, but contain a significant amount of error. Refining these estimates 
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has two parts: mosaicing and registration [Coorg et al., 1998, Coorg, 1998]. 
Mosaicing exploits the constraint that all of the images in a given node share 
the same center of projection and form (somewhat more than) a hemisphere. 
The relative orientations of a node are optimized by correlating the overlap- 
ping regions of adjacent images. A sample of the spherical mosaic produced is 
shown in Figure 1-7. Once mosaicing is completed a node can be treated as 
a single image with a single position and orientation. Registration identifies 
a few prominent features per node and performs a global optimization which 
refines the location and orientation of each node. 



1.3 Thesis Overview 

Figure 1-8 shows an overview of the reconstruction pipeline described in this 
thesis. The calibrated i magery descri bed in the last section serves as input to 
the pipeline. The left hand column shows the major steps of our approach; the 
right hand side shows example output at various stages. The output of the 
pipeline is a textured three-dimensional model. Our approach can be generally 
characterized as hypothesize and test. We hypothesize a surfel and then test 
whether it is consistent with the data. Chapter 2 presents our basic approach. 
To si mpl ify the presentati on i n thi s chapter we assume perfect data (i .e. perfect 
calibration, constant illumination, diffuse surfaces, and no occlusion or image 
noise) and use a synthetic dataset to demonstrate the algorithm. Chapter 3 
extends the theory to cover camera calibration error, variations in illumina- 
tion, occlusion, and image noise. Chapter 4 explores detecting and localizing 
surfels. We demonstrate the ability of our algorithm to compensate for cam- 
era calibration error, variations in illumination, occlusion, and image noise. 
We examine its ability to detect and localize three-dimensional surfaces from 
a significant distance (both in position and orientation) using a purely local 
algorithm. Noisy data causes sometrue positives to be missed and compensat- 
ing for noisy data increases the number of false positives recovered. Chapter 5 
introduces several techniques to remove false positives and fill in the missing 
parts. Camera updates and grouping surfels into surface impose global con- 
straints which remove nearly all of the false positives. Growing surfaces fills 
in most of the reconstruction. Finally we discuss fitting simple models and ex- 
tracting textures. We present the algorithms as well as the results of applying 
them to a large dataset. 

1.3.1 Definitions 

Some of the more general or fundamental terms used throughout this thesis 
are defined below. Others will be defined as needed. 



1.3. THESIS OVERVIEW 



27 




Detection 

and 

Localization 



I 



Camera Calibration 

Updates and 

Outlier Rejection 



Grouping Surfels 



Growing Surfaces 



Building Models 

and 

Extracting Textures 









dtti-nV, 


r*Jj 


i ^^f^ m 




• W^^r ■ 






Figure 1-8: Overview of Thesis. 
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• Internal Parameters: The camera parameters which govern the image 
formation process. The exact parameters depend on the camera model 
used, but usually include: focal length, aspect ratio, and principle point. 
See Appendix A for a more complete description. 

• External Parameters: The camera parameters which govern what, in 
an absolute coordinate system, is imaged - the location of the center of 
projection and the orientation of the optical axis. See Appendix A for a 
more compl ete descr i pt i on . 

• Calibrated Image: An image for which both the internal and external 
calibration are known. Also referred to as a pose image. 

• Node: A set of images acquired from the same location. I n other words, 
all of the images in a node share the same center of projection. The 
images of a node typically have orientations which tile a hemisphere or 
moretCoorg et a I., 1998]. 

• Surfel: A small planar surface patch or surface element. This definition 
differs from Szeliski's [Szeliski and Tonnesen, 1992] in that it refers to a 
finitesized patch which includes both geometry and texture. 

• Units: All of the datasets used in this thesis share a common global coor- 
dinate system in which 1 unit equals 1/10 foot. 

The pin-hole camera model is used throughout this thesis. As noted in Ap- 
pendix A it is a linear model and is not capable of modeling nonlinear dis- 
tortion. Images which have large amounts of distortion can be unwarped as 
a preprocessing step to remove the nonlinear portions, however this was not 
necessary for the datasets used in this thesis. 



C hapter 2 

The Basic Approach 



In this chapter we describe our basic approach. To simplify the presentation 
we assume perfect data (i.e. perfect calibration, constant illumination, diffuse 
surfaces, and no occlusion or image noise). These assumptions are relaxed in 
the following chapter. The approach presented here was initially inspired by 
the work of Bolleset al. [1987]. The notation used in this chapter is defined in 
Tabl e 2. 1 and Secti on 2.1 revi ews epi pol ar geometry and epi pol ar-pl ane i mages. 
Wethen definea construct called an epi polar image (Secti on 2.2) and show how 
it can be used to reconstruct three-dimensional information from large sets of 
images (Section 2.3). LikeTsai [1983] and Okutomi and Kanade[1993] we de- 
fine a cost function that is applied across multiple images, however, we do not 
evaluate the cost function in image space. Instead, like Collins [1996], Seitz 
and Dyer [1997], and Kutulakos and Seitz [1998b, 1998a] we perform match- 
ing in three-dimensional space. We also model occlusion and the imaging pro- 
cess similar to Cox [1996] and Belhumeur [Belhumeur and Mumford, 1992, 
Belhumeur, 1993, Belhumeur, 1996]. There are several important differences, 
however. An epi polar image is valid for arbitrary camera positions and over- 
comes some forms of occlusion. Where three-dimensional information cannot 
be recovered using purely local information, the evidence from the epi polar 
image provides a principled distribution for use in a maximum-likelihood ap- 
proach (Section 2.4) [Duda and Hart, 1973]. Finally, we present the results of 
using epi polar images to analyze a set of synthetic images (Section 2.5). 



2.1 Epi polar Geometry 

Epi polar geometry provides a powerful constraint for identifying correspon- 
dences. Given two cameras with known centers Ci and C 2 and a point P in the 
world, the epi polar plane n c is defined as shown in Figure 2-1. P projects to 
p 1 and p 2 on image planes U} and n 2 respectively. The projection of II, ontollj 1 
and nf produce epi polar lines^ and l\. This is the essence of the epi polar con- 
straint. Given any point p on epi polar I i ne 4 i R imager^ 1 , if the corresponding 
point is visible in image n 2 , then it must Neon the epi polar line^. 

29 
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Pj Absol ute coordi nates of the j th surfel . 

rij Orientation of the j th surfel. 

Cj Center of projection for the i th camera. 

n? I mage plane for the i th camera. 

r Cal i brated i mage, r = (nj, Q) . 

I Set of calibrated images, {f}. 

p* I mage poi nt. Projection of P, onto n?. 

n^ The /c th epi pol ar pi ana 

e*>* E pi pol ar line. P roj ecti on of n^ onto n{ . 

p* Base image point. Any point in any image. 

Pj Set of proj ecti ons of P., , |pj } . 

p^ Set of projections of Pj for which the projection is in the 

front half space, {p* C~P* • n< o}. 

C* Base camera center. Camera center associated with p*. 

C Set of camera centers, {CJ. M ay or may not i ncl ude C* 

n* Base image. Contains p*. 

Ili Set of image planes, {nj}. May or may not include n* 

t Baseline. 3D line passing through p* and C*. 

t/ Epi polar line. Projection of t ontoIP ; . 

t e Set of epi polar lines, {t/}. 

SP k Epi polar plane image. Constructed using n^. 

S x Epi polar image constructed using p*. x indexes all p*'s. 

T(P,I) Function which projects P onto I, eg. T(P^P) = pj. 

V(P,n,I) Ideal i mage values (no noise or occlusion) at T(P, I). 

F{x) Function of the image at points (eg. pixel values, 

features, etc). 

X(x 1 ,x 2 ) Matching function. Quality of match between x 1 andx 2 . 

v{j,a) Match quality. Used to analyzed. Alsoz/(P i ,n i ) and v(Sj] 

d(p*) Depth at image point p\ Distance from d to the closest 

— ► 

actual surface in the direction Cip J . If low confidence 
or unknown, then oo. 
o(p*) Orientation at image point p*. Orientation of to the 

closest actual surface in the direction Qp\ 
G(/i, a 2 ,x) Gaussian distribution with mean ji and variance a 2 

evaluated at x. 
{E\C} Set of al I E's such that C is true. 
p(P | C) Probability of P conditioned on C. 
P1P2 Unit vector in the direction from Pi toP 2 . 

Mi Modeled object. Previously reconstructed. 



Table 2.1: Notation used for the basic approach. 
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Figure 2-1: E pi polar geometry. 



Bolles et al. [Bolles et al., 1987] used the epi polar constraint to construct 
a special image which they called an epi polar-plane image. As noted earlier, 
an epi polar line 4 contains all of the information about the epi polar plane H 
that is present in thei th image Uf. An epi polar-plane image is built using all of 
the epi polar lines^ (the boldface symbol denotes a set, i.e. U^}) from a set 
of images 11; which correspond to a particular epi polar plane 1^ (Figure 2-2). 
Since all ofthelines^ in an epi polar-plane image £P k are projections of the 
same epi polar plane n^, for any given pointp in SP k , if the corresponding point 
in any other image nj is visible, then it will also be included in £P k . Bolles, 
Baker and Marimont exploited this property to solve the correspondence prob- 
lem for several special cases of camera motion. For example, with images taken 
at equally spaced points along a linear path perpendicular to the optical axes, 
corresponding points form lines in the epi polar-plane image; therefore finding 
correspondences reduces to finding lines in the epi polar-plane image. 

For a given epi polar plane nj;, only those images whose camera centers lie 
on n^ (JQ Cj-n|: = o}) can be included in epi polar-plane imaged. For ex- 
ample, using a set of images whose camera centers Neon a plane, an epi polar- 
plane image can only be constructed for the epi polar plane containing the cam- 
era centers. In other words, only a single epi polar line from each image can be 
analyzed using an epi polar-plane image. In order to analyze all of the points in 
a set of images using epi polar-plane images, all of the camera centers must be 
collinear. This can be a serious limitation. 



2.2 E pi polar I mages 

For our analysis we define an epi polar imaged which is a function of a set of 
images and a point in one of those images. An epi polar image is similar to 
an epi polar-plane image, but has one critical difference that ensures it can be 
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Figure 2-2: E pi polar-plane i mage geometry. 



c. 




Figure 2-3: E pi polar i mage geometry. 
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constructed for every pixel in an arbitrary set of images. Rather than use pro- 
jections of a single epi polar plane, we construct the epi polar image from the 
pencil of epi polar planes defined by the linef through one of the camera cen- 
ters C* and one of the pixels p* in that imagen* (Figure 2-3). n* is the epi polar 
plane formed by t and thei th camera center Q. Epi polar line 4 contains all 
of the information about t present in n?. An epi polar-plane image is composed 
of projections of a plane; an epi polar image is composed of projections of a line. 
The cost of guaranteeing an epi polar image can be constructed for every pixel is 
that correspondence information is accumulated for only one point 1 p*, instead 
of for an enti re epi polar I i ne. 




Figure 2-4: Set of points which form a possible correspondence. 



To simplify the analysis of an epi polar image we can group points from the 
epi polar lines according to possible correspondences (Figure 2-4). Pi projects 
top 1 ! in nj; therefore Pi ({pi}) has all of the information contained inH about 
Pi. There is also a distinct set of points P2 for P 2 ; therefore p., contains all of 
the possible correspondences for Pj. If P., is a point on the surface of a physical 
object and it is visible in both n? and n*, then measurements taken atp* should 
(under thesimpleassumptions of this chapter) match those taken atp* (Figure 
2-5). Conversely, if P, is not a point on the surface of a physical object then the 
measurements taken atp* are unlikely to match those taken atp* (Figures 2-6 
and 2-7). Epi polar images can be viewed as tools for accumulating evidence 
about the possible correspondences of p*. A simple fundi on of j is used to build 



Y lo simplify the presentation in this chapter the discussion will focus on points, however 
(oriented) points and surfels are interchangeable. 
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an ordered set so that {P.,} is a set of points along t at increasing depths from 
the image plane. 





Figure 2-5: P, is a point Figure 2-6: Occlusion be- Figure 2-7: Inconsistent 



on a physical object. 



tween C,- and P.- 



background 



2.3 Basic Algorithm 

An epi polar imaged is constructed by organizing 

{^(pp \F(-) is a function of the image} 

into a two-dimensional array with i and j as the vertical and horizontal axes 
respectively. Rows in £ are epi polar lines from different images; columns form 
sets of possible correspondences ordered by depth 2 (Figure 2-8). The quality 
v(j) of the match 3 between column j and p* can bethought of as evidence that 
p* is the projection of P., and j is its depth. Specifically: 



v{ J ) = Y J X{Hv l j ),Hv*)), 



(2.1) 



where F{-) is a function of the image and X{-) is a function which measures 
how well Fivtj) matches ^"(p*). A simple case is 



F{x) = i ntensity val ues at x 



and 



X(xi,x 2 ) 



\xi — x 2 \ 



2 The depth of Pj or distance from C* to Pj can be trivially calculated from j, therefore we 
consider j and depth to be interchangeable. We further consider depth and three-dimensional 
position to be equivalent since we use calibrated images and the three-dimensional position 
can be trivially calculated from p* and the depth. 

interestingly, v(j) is related to the tomographic algorithm [Ramachandran and Lakshmi- 
narayanan, 1971, Geringand Wells, 1999]. 
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S(I,l*,p* 



i (image) 



Figure 2-8: Constructing an epi polar image. 



Real cameras are finite, and pj may not be contai ned i n the i mage nj (p) g n? 
Only terms for which p* e n? should be included in (2.1). To correct for this, 
v(j) is normalized, giving: 



"(J) 



y: ^(^(pj),^(p*)) 

i Ipjsn? 

i|p}en? 



(2.2) 



Ideally, i/(j) will have a sharp, distinct peak at the correct depth, so that 



argmax(z/(j')) 

3 



the correct depth at p* 



As the number of elements in p^ increases, the likelihood increases that v{j) 
will be large when P^ lies on a physical surface and small when it does not. 
Occlusions do not produce peaks at incorrect depths or false positives 4 . They 
can however, cause false negatives or the absence of a peak at the correct depth 
(Figure 2-9). A false negative is essentially a lack of evidence about the correct 
depth. Occlusions can reduce the height of a peak, but a dearth of concurring 



^Except possibly in adversarial settings. 
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images is required to eliminate the peak. Globally this produces holes in the 
data. While less then ideal, this is not a major issue and can be addressed in 
two ways: removing the contribution of occluded views, and adding unoccluded 
views by acquiring more images. 



Mi 



M 2 




Figure 2-9: False negative caused by occlusion. 




Figure 2-10: Exclusion region (grey) for surfel located at P, with normal n,- 



A large class of occluded views can be eliminated quite simply. Figure 2-10 
shows a surfel located at point P, with normal n,-. I mages with camera centers 
in the hashed half space cannot possibly have viewed P,. n,- is not known a 
priori, but thefact that P., is visible in nf limits its possible values. This range 
of values can then be sampled and used to eliminate occluded views from v(j). 
Let a bean estimate of the normal 5 n, and QP_) betheunit vector along the ray 



5 Actually, a is the azimuth and the elevation is assumed to be 0. This simplifies the presen- 
tation without loss of generality. 
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fromCj toPj,thenPj can only be visible in nj if CjPj-a < 0. We denote the set of 
image points with meet this visibility constraint asp". At a fundamental level 
an estimated normal is associated with each reconstructed three-dimensional 
point (or surfel). For single pixels, the estimated normal is very coarse, but its 
existence is very useful for grouping individual points (or surfel s) into surfaces. 
Neighborhood information, introduced in Chapter 3 as small planar patches, 
greatly improves the accuracy of the estimated normals. 

If the volume imaged by H is modeled (perhaps incompletely) by previous 
reconstructions, then this information can be used to improve the current re- 
construction. Views for which the depth 6 at pj, or d(pj), is less than the distance 
from d to Pj can also be eliminated. For example, if Mi and M 2 have already 
been reconstructed, then the contributions of nj where i e {1,2,3} can be elim- 
inated from v(j) (Figure 2-9). In addition, we weight the contributions based 
on their forshortening. The updated function becomes: 



v(j, a) = 


E(^-«)*(-F(p5).^(p*)) 

ieQ 

ieQ 


where 

Q=< 


f 

i 


d( P }) > HC-P.-H 2 J 


Then, if sufficient evidence exists, 


argmax(V(j, 

3, a 


a) 


f j = depth at p* 
1 a = esti mate of n,- 



(2.3) 



2.4 Probabilistic Formulation 

Probabilities can be used to formulate the notion that v(j) and v(j,a) should 
be large when P 7 lies on a physical surface and small otherwise. Given image 
data I, q the probability that P., lies on a physical surface, has depth j and 
orientation a, and is visible in n* can be written formally as 

g = p(d(p*)=j,o(p*) = a|I). (2.4) 

Using Bayes' rule gives 

p(I | d(p*) = j, o(p*) = a)p(d(p*) = j, o(p*) = a) 



q 



P(I) 



6 Distance from Q to the closest previously reconstructed object or point in the direction 
Qp}. 
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Of all the image points in I only p^ depend upon a and j. The rest can be 
factored out yielding 

PiPj I d(p*) = j, o(p*) = a )p(d(p*) = j, o(p*) = a) 

1 = ?~^ (Z5) 

P\Pj) 

log(g) = log(p(p? | d(p*) = j, o(p*) = a)) + (2.6) 

Iog(p(d(p*) = j,o(p*) = a)) - log(p(p?)) 

If we assume that the measured pixel values ^(pj) contain Gaussian noise 
with a mean value of .F(p*) and a variance of a 2 and that individual pixel mea- 
surements are independent then the first term becomes 

~\ E {(HV)) ~ W)) V + log (2vra 2 )) (2.7) 

i 

which is very similar to v(j, a) (Equation 2.3). The next term of Equation 2.6 
is a prior on the distribution of depths and orientations. If a partial model 
already exists, it can be used to estimate these distributions. Otherwise we 
make the standard assumption that all j'sand a's are equi -probable. 

The last term is more challenging; how doweestimatep(p")? Applying the 
same i ndependence assu mpti on as above yi el ds 

p(p?) = £log(p(pj)). 

i 

We could assu me that all values of T(-) areequi-probablein all images, making 
the denominator irrelevant, and end up with a matching function equivalent 
to Equation 2.3. Another approach is to use nf toestimatep(pj). p(p*-) can be 
thought of as a significance term. The significance of a match between ^"(p*) 
and .F(p*) is inversely proportional top(pj-). We use a nearby (in both position 
and orientation) image nf instead of nf toestimatep(pj). pj, C u and C k define 
an epipolar plane nf. To estimate the distribution we consider epipolar line 
e** (projection of n** onto image plane nf). e** contains all of the possible 
matches for p*. in nf and we use this set of possible matches and a mixture of 
Gaussians model toestimatep(p}) giving: 

£ G(^(p^),a 2 ^( P )) 



p(p}) = ^ =- . (2.8) 



P6^ !,t 



Since p(pj) depends only upon pj and nf, it can be computed once as a pre- 
processing step. Figure 2-11 shows two probability images produced in this 
manner. Black represents low probability and white high. p(pj) can bethought 
of as a uniqueness term -p* matching p* is only significant ifp(p^) is low. 
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Figure 2-11: Probability images for synthetic data. 

Our approach to estimating p(d(p*) = j, o(p*) =a |I) is conservative. It does 
not consider subsets of p" or allow individual elements to be tagged as occluded 
or outliers. Nevertheless valuable insight was gained fromexaminingp(p^). In 
the next section results both with and without considering^^) are presented. 

2.5 Results 




-i. ti i-i 



Figure 2-12: Example renderings of the model. 
Syntheti c imagery was used to explore the characteristics of u(j) and u(j, a) 
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A CAD model of Technology Square, thefour-building complex housing our lab- 
oratory, was built by hand. The locations and geometries of the buildings were 
determined using traditional survey techniques. Photographs of the buildings 
were used to extract texture maps which were matched with the survey data. 
This three-dimensional model was then rendered from 100 positions along a 
"walk around the block" (Figure 2-12). From this set of images, a If and p* 
were chosen and an epi polar imaged constructed. £ was then analyzed using 
two function, Equations 2.2 and 2.3where 



and 



T(x) = rgb(x) = [r(x),g(x),b(x)] (2.9) 



#([ri, si, &i],[r 2 , 02,62])= (2.10) 



(n - r 2 ) (01 - 02) (61 - 62) 



2- 



°r °l °l 



r(x), g(x), and b(x) are the red, green, and blue pixel values at point x. of, a\, 
and <j\ are the variances in the red, green, and blue channels respectively. For 
the synthetic imagery used in this chapter we set of = a\ = a\ = 1. Elsewhere, 
the variances are estimated during the internal camera calibration and are 
assumed to remain constant. 

Figures 2-13 and 2-14 show a base image If with p* marked by a cross. 
Under n* is the epi polar imaged generated using the remaining 99 images. 
Below £ is the matching function v(j) (Equation 2.2) and v{j,a) (Equation 2.3). 
The horizontal scale, j or depth, is the same for £ , u(j) and v(j,a). The vertical 
axis of £ is the image index, and of v(j, a) is a coarse estimate of the orientation 
a at Pj. The vertical axis of ^(j) has no significance; it is a single row that has 
been replicated to increase visibility. To the right, v(j) and v(j,a) are also 
shown as two-dimensional plots 7 with the correct depth 8 marked by a line. 

Figure 2-13a shows the epi polar image that results when the upper left- 
hand corner of the foreground building is chosen as p*. Near the bottom of £, 
t c is close to horizontal, and p* is the projection of blue sky everywhere except 
at the building corner. The corner points show up in £ near the right side as 
a vertical streak. This is as expected since the construction of £ places the 
projections of P., in the same column. Near the middle of £ , the long side to 
side streaks result because P^ is occluded, and near the top the large black 
region is produced because p* £ n?. Both v(j) and u(j,a) have a sharp peak 9 
that corresponds to the vertical stack of corner points. This peak occurs at a 



Actually, £ tt u(j, a)/£ Q 1 is plotted for v {j, a). 
determined by intersecting £* with the model. 
9 White indicates the best match, black the worst. 
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Figure 2-13: U*,p*,8,u(j) and i/(j, a) (Part I). 
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Figure 2-14: n?, p*, E, v(j) andi/(j,a) (Part II). 
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depth of 2375 units 10 0' = 321) for v{j) and a depth of 2385 (j = 322) for u(j, a). 
The actual distance to the corner is 2387.4 units. The reconstructed world 
coordinates of p* are [-1441,-3084,1830] and [-1438,-3077,1837] respectively. 
The actual coordinates 11 are [-1446, -3078, 1846]. 

Figure 2-13b shows the epi polar image that results when a point just on 
the dark side of the front left edge of the building is chosen asp*. Again both 
v(j) and v{j,a) have a single peak that agrees well with the actual depth. This 
time, however, the peaks are asymmetric and have much broader tails. This is 
caused by the high contrast between the bright and dark faces of the building 
and the lack of contrast within the dark face. The peak in v(j,a) is slightly 
better than the one in v(j). 

Figure2-13c shows the epi polar image that results when a point just on the 
bright side of the front left edge of the building is chosen asp*. This time v(j) 
and v(j,a) are substantially different. v{j) no longer has a single peak. The 
largest peak occurs at j = 370 and the next largest at j = 297. The peak at 
j = 297 agrees with the actual depth. The peak at j = 370 corresponds to the 
point where t exits the back side of the building. Remember that v(j) does 
not impose the visibility constraint shown in Figure 2-10. u(j, a), on the other 
hand, still has a single peak, clearly indicating the usefulness of estimating a. 

In Figure 2- 14a, p* is a point from the interior of a building face. There 
is a clear peak in v(j,a) that agrees well with the actual depth and is better 
than that in u(j). In Figure 2-14b, p* is a point on a building face that is 
occluded (Figure2-9) in a number of views. Both v(j) and v(j,a) produce fairly 
good peaks that agree with the actual depth. In Figure 2- 14c, p* is a point 
on a building face with very low contrast. In this case, neither v(j) nor v(j,a) 
provide clear evidence about the correct depth. The actual depth occurs at 
j = 387.5. Both v(j) and v(j,a) lack sharp peaks in large regions with little or 
no contrast or excessive occlusion. Choosing p* as a sky or ground pixel will 
produce a nearly constant v(j) or u(j, a). 

Figures 2-15 and 2-16 show the result of running the algorithm on synthetic 
data. A grey circle and a black line are used mark the location and orientation 
of each image in the dataset. There and y coordinates have been divided by 
1000 to simplify the plots. For each pixel in the set of images shown in Figure 
2-12 an epi polar image was constructed and the point with maximum v(j,a) 
recorded. Points for which u(j, a) is at least -2 and |p,| is at least 5 are plotted. 
Global constraints such as ordering or smoothness were not imposed, and no 
attempt was made to remove low confidence depths or otherwise post-process 
the maximum of u(j, a). In Figure 2-15 the grey levels correspond to the density 
of reconstructed points, with black the most dense. In Figure 2-16 the grey 



10 lunit is 1/10 of afoot. 

n Some of the difference may be due to the fact that p* was chosen by hand and might not be 
the exact projection of the corner. 
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Figure 2-15: Density of reconstructed points. 
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Figure 2-16: Orientation of reconstructed points. 
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Figure 2-17: Distribution of errors for reconstruction. 
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Figure 2-18: Two views of the reconstruction. 
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levels correspond to orientation. The grey scale axis to the right is in units of 
7r/2. The outlines of the four buildings are clearly visible and the orientations 
are roughly correct. 

What is difficult to appreciate from these plots is the relative numbers of 
points contributing to the black lines compared to the grey outliers scattered 
about. Figure 2-17 shows the error distribution for the reconstruction. Plotted 
val ues i ndicate the percentage of reconstructed poi nts whi ch have no morethan 
the indicated error. The error measure is 

||Pj (reconstructed) — P.,- (actual) || / ||C* — Pj (actual) || . 

Notethat > 85% of the reconstructed points are within 1% of thecorrect value. 
Figure 2-18 shows the lower left corner of the bottom building as viewed from 
[-1600, -4000, 1100] and [-2000, -3500, 1200]. The results are rendered as ori- 
ented rectangular surfaces using the reconstructed position and estimated ori- 
entation. The size is set so that the projection in n* is 1 pixel. For clarity, 
the reconstructed points have been down-sampled by 3 in each direction. The 
recovered orientations may not be good enough to directly estimate the under- 
lying surface, but they are good enough to distinguish between surfaces. This 
idea will play an important role in Chapter 5. 

Next we explore several modifications to the basic algorithm described in 
Section 2.3. The synthetic images in Figure 2-12 were rendered with direc- 
tional lighting in addition to global lighting. As a result, image brightness 
varies with viewing direction. To compensate we add a brightness correction 
term 7 to the matching function 12 : 



X{[ri,gi,b 1 \,[r2,g2,b 2 \) 



7 



/ (TTi -r 2 ) 2 (Tgi - g 2 f (761 -b 2 )^ 
\ <*\ °l o\ 

rir2 1 9192 1 fei&2 



l>1 



- 2 a 2 

1 I i>l I "1 

~T ~r — 2 ~r ~% 

'r g b 



We also consider two variations of the probabilistic formulation developed in 
Section 2.4. In one we limit log(p(p")) > -15 and the other we do not. The 
upper left image shown in Figure 2-12 was selected as H*. 39% or 67,383 pix- 
els are reconstructive (i.e. non-sky and non-ground). Each pixel in If was 
reconstructed using the six variants shown in Table 2.2. Variants which in- 
clude brightness compensation are annotated with "w/ be" and those without 
"w/o be". Variants which use the probabilistic formulation are annotated "+P" 
and thosethat limit the log probability have "> =-15" added. Figure 2-19 shows 
error distributions for the variants. Interestingly, variant 6 has the best error 



12 This is similar to normalized correlation, however we impose nonlinear constraints on 7. 
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Figure 2-20: Error distribution after adding noise. 
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distribution, but misses many reconstruct ble points. Conversely, variant 1 re- 
constructs nearly all of the points, but sacrifices some accuracy. Variant 3 is 
the overall best. 



Variant 


Label 


Points Reconstructed 


Number 


Percent 


1 


M w/bc 


65,582 


97 


2 


M -+P w/ be 


63,579 


94 


3 


M-rf D >=15 w/bc 


62,564 


93 


4 


M w/obc 


34,710 


52 


5 


M -+P w/o be 


24,281 


36 


6 


M-rf D >= 15 w/obc 


23,485 


34 



Table 2.2: Performance of algorithm variants. 

Finally, we consider the effects of camera calibration noise. Ultimately it is 
hoped that Argus will achieve positional accuracy on the order of a few centime- 
ters and pose accuracy on the order of a few mi Mi radians. Argus has not yet 
achieved this design goal, but we use it as a reference. Uniformly distributed 
noise was added to each coordinate of the position and to the orientation in 
increments of 2" and 0.1° respectively. The results are shown in Figure 2-20. 

2.6 Discussion 



This chapter defines a construct called an epi polar image and then uses it to 
analyze large sets of synthetic data. This analysis produces an evidence versus 
position and surface normal distribution that in many cases contains a clear 
and distinct global maximum. The location of this peak determines the po- 
sition. The algorithm presented uses only local calculations and lends itself 
nicely to parallelization. 



Chapter 3 

The Challenge of Noisy Data 



This chapter extends the approach presented in the previous one to work with 
noisy data. Specifically data which contains 1) camera calibration error; 2) 
significant variations in illumination; and, 3) complex occlusion (eg. building 
viewed through tree branches) and image noise. With these modifications, our 
approach shifts from looking for sets of matching points to searching for sets 
of highly correlated regions. The fact that at a fundemental level we match 
surfels (small planar regions with position, orientation, and texture) is made 
explicit in this chapter. Theaddition of neighborhood information produces sev- 
eral benefits, such as accurate normal estimation, but also makes the epi polar 
i mage construct £, described in Chapter 2, less useful. £ can no longer be used 
to directly visualize possible correspondences. Hypothesizing points along t 
also becomes less attractive. This point wil be addressed furter in Chapter 4. 
Some additional notation used in this chapter is defined in Table 3.1. 



3.1 Camera Calibration Errors 

The epi polar constraint exploited in Chapter 2 relies on accurate camera cali- 
bration. If the calibration is not accurate then the constraint that correspond- 
ing points must lie on a pair of epi polar lines (eg. 4 and t\ in Figure 2-1) no 
longer holds. If the calibration error can be bounded then the epi polar lines 
become epipolar stripes. Consider perturbing Q and C 2 in Figure 3-1 about 
their centers of projection; the intersection of lie with Uj and nf sweeps out 
l\ and l\. The area of the epipolar stripe is proportional to the calibration 
error. What was a one-dimensional search with perfect calibration is now a 
two-dimensional search. As shown in Figure 3-2, the actual projection of P is 
not p. If the bounded camera calibration error 5 is known thenf> can be calcu- 
lated andp ep. The effect can be modeled as a displacement [u,v] of the actual 
projection. 

Calibration error complicates the epipolar image £ and its analysis de- 
scribed in Sections 2.2 and 2.3. p^ does not contain all of the information in 
ft? about P^; we must consider p^. I nstead of an array (indexed by j and a) of 
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Figure 3-1: E pi polar stripes. 





Figure 3-2: Effect of camera calibration error. 
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X"] 



Sj j th surfel. 

y x Sj Point in Sj. x and y are indices intoS,. gS^ = Pj. 

sj P roj ect i on of Sj onto ft? . 

%j Set of proj ecti ons of Sj , |s* \ . 

Point in sj. a; and y are indices intosj. [jsj = pj. 

s* Basesurfel. {^s*}. May contain estimated values. 

8j M atch i n g set of su rf el s. [ sj sj matches s* } . 

Q Estimated center of projection for thei th camera, 

ft? Estimated image plane for thei th camera. 

r Cal i brated i mage (esti mated cal i brati on), r = (nj, Q) . 

pj I mage poi nt. P roj ecti on of P., onto nj . 

pj Set of proj ecti ons of P.,- . {pj } ■ 

pj Noisy Projection of P., onto ft*. Area that 

Pj projects to given bounded camera calibration error. 

£pj Point inpj. u and v are indices intopj. 

Pj %v l j where Ui and ^ produce the best match. 

pj Set of matching image points. {pj pj matches p*}. 

Pj Reconstructed world point. Estimated using p.,. 

Pj Set of noi sy proj ecti ons of Pj , {pj } . 

n£ The k th esti mated epi pol ar pi ane. 

l k / Epi polar stripe. Noisy projection of P^ onto ft?. 

S i Bounded camera calibration error for f 

T(P,I,«J) Noisy projection function, eg. T(Pj,P,^) = pj. 

NoteT(P i ,?,<J i ) = pj. 



Table 3.1: Notation used for handling noisy data. 

individual pixels, £ is now composed of two-dimensional regions. Equations 2.2 
and 2.3 become 

KJ) = l2 ~ y~ x (3-D 



i pj en? 



and 



^C.P, • a] max A^Gp}),^)) 

u(j, a) = ** =^ (3.2) 

£QP, " a 

«eS 
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where 

Q = 



p] e uf 






V) G P? 




d(pj) > 


Q- 


-Pj 



Equation 2.8 must also be modified. E pi polar stripelf fe is used to estimate the 
distribution instead of epi polar lineup producing 



. Tffc ?'. . A' 



p(pj) = ^ =- . (3.3) 

I n the last chapter, if a match was found the reconstructed point was auto- 
matically Pj. This is no longer the case. p\, represents our best estimate for a 
set of corresponding points and in and v { are shifts which correct for the cali- 
bration error. What world point P, gave risetop^? It almost certainly wasn't 
Pj. Assuming that the calibration error can be modeled as zero mean additive 
Gaussian noise then the best estimate of R, is the one which minimizes the 
sum of squared calibration errors. Two possible measures of the calibration 
error are the distance (in three-dimensional space) between P, and the line of 
sight through pj, 



argmin V Qp} x ((P, - Q) x Qpj) • (P, - Q) (3.4) 

J pjePi 

and the distance (in two-dimensional space) between pj and the projection of 

P , T(P-,P) 

argmin V |p» - T(P,-,P)| 2 - (3-5) 

p . — • J 
3 pjePi 

The shifts {[^i,^]} introduce additional degrees of freedom. One benefit of this 
is that Pj can be reconstructed from more than one P.,, allowing P, to be recov- 
ered even if P., is not very close. This property will be exploited in Chapter 4. 
The additional degrees of freedom also admit additional false positives. Figure 
3-3 shows one example which occurs because similar looking points (the upper 
left corner of a black square) lie within each j^ (the grey circles). One way 
to limit the number of false positives is to constrain the shifts. All of the shifts 
applied to points in a given image should be consistent with a single camera 
correction. This will be explored further in Chapter 5. 

3.2 Variations in Illumination and Reflectance 

The matching function proposed in Equation 2.10 assumes controlled lighting 
and diffuse reflectance. Only the overall brightness is allowed to change. Nat- 
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Figure 3-3: False positive induced by compensating for camera calibration er- 
ror with shifts. 



ural lighting can vary significantly (eg. time of day, season, weather, etc.) not 
only in brightness, but also chromaticity Real scenes are also not restricted to 
objects which reflect diffusely. Both of these make it difficult to compare im- 
ages of the same region taken under differing illumination conditions and from 
different viewpoints. 

A simple linear model can be used to correct images taken under different 
viewing conditions (illumination and/or reflectance). The radiance R arriving 
at an image sensor from a given region is simply the irradiance/ at the region 
multiplied by its reflectance A;: 

R = kl. (3.6) 

and the value measured by the image sensor is: 



f(R) + d. 



(3.7) 



f(-), an invertible function, and b, an offset, are functions of the image sen- 
sor 1 . Figure 3-4 shows a region imaged under different conditions. Clearly, if 
hh/hh is constant for some region then R 2 /R\ will also be constant over the 
same region. From Equations 3.6 and 3.7 it follows that if for some region 



Ioko = Lhic 



(3.8) 



Equation 3.7 i s commonly called the image sensor's intensity response function [Voraetal., 
1997a, Vora eta I., 1997b]. 
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Figure 3-4: Region imaged under different conditions. 



and 

then 

or 



f(xy) = f(x)f(y) 



v 2 = f(c)vi - f(c)di + d 2 



v 2 = mv i + d. 



(3.9) 



(3.10) 



Equation 3.9 is true for typical CCD cameras. Equation 3.10 applies to both 
diffuse and non-diffuse surfaces and can be calculated independently for each 
color channel. We often refer torn and d as an illumination correction. It 
allows us to remap image values obtained under one condition to another and 
then match them. Substituting Equation 3.10 into Equation 2.10 gives: 



X{[ri,gi,b 1 \,[r2,g2,b 2 \) = 

( '((m r ri + dr) - r 2 f {{m g gx + d s ) - g 2 ) 2 {{m h b x + d s ) - b 2 \ 



(3.11) 



at 



at 



at 



The condition given by Equation 3.8 requires that the product of irradiance 
and reflectance vary similarly for changes in viewing conditions at each point 
x in the region 11: 

Wxen:I 2 r 2 /I iri =c. (3.12) 

In practical terms, this means that TZ must be planar (or well approximated 
by a plane), uniformly lit, and contain materials with similar reflectance func- 
tions. If we are willing to assume that only a limited number of materials 
appear in each region, then using clustering techniques to calculate the cor- 
rection removes the requirement that the reflectance functions be similar. For 
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example, in the next chapter, a region containing a window with specular re- 
flection and concrete is successfully corrected with this technique. The spectral 
content of theirradiance may be changed arbitrarily from one condition to an- 
other and the region may contain different colored materials. If Equation 3.9 
is not valid and the intensity response fund ion is known, then image values 
can be converted to radiance and the remapping done there. 

A single pixel was sufficient to calculate v(-) as presented in Sections 2.3 
and 3.1. This is no longer the case once we consider illumination. In order to 
calculate m and b at least 2 pixels are needed and for good results more should 
be used. Consider matching images of a surfel (small planar region) S, with 
orientation n, located at P.,. Since orientation is intrinsic to a surfel, Equation 
3.1 is no longer useful. We can rewrite Equation 3.2 as 

W^-n,)max £ X (Hill*)) , H^)) I £ 1 
u(Sj) = v J : L__ J J ' . (3.13) 

£QP, . n . 

ieQ 

In order toobtain a match, S.,- must beclose(both in position and orientation) to 
a physical surface and Equation 3.12 must be satisfied. Correcting for viewing 
conditions allows us to match images taken under different illumination con- 
ditions and from different viewpoints, but also admits false positives. Figure 
3-5 shows one example. One way to limit the number of false positives is to 
constrain the correction (mr,d r ,m g ,d g ,m h ,d b ). This is explored in Section 4.2. 

There are several drawbacks to matching surfels instead of individual pix- 
els: the cost of computing the match is higher; the shifts introduced in Section 
3.1 vary with image position and the distance to the imaged point and there- 
fore are not constant throughout the surfel; and, calculating p(pj-) is problem- 
atic. For small calibration errors and small regions it is reasonable to assume 
that the shifts are piecewise constant and because of the additional informa- 
tion content of a surfel, p(pj) is less important. On the positive side, the image 
data of a surfel can be used to efficiently calculate in and v { and accurately 
estimate its normal (Sections 4.2 and 4.3). 

3.3 Occlusion and I mage Noise 

Occlusion poses another difficult problem. Figure 3-6 shows two examples typ- 
ical of urban environments. On the left, the foreground building completely 
occludes two columns of windows. This view provides no information about 
the occluded region and if enough other views are completely occluded we will 
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Figure 3-5: False positive induced by correcting for viewing conditions. 
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Figure 3-6: Hard (left) and soft (right) occlusion. 
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not be able to reconstruct that portion of the background building. The exam- 
ple on the right is potentially more difficult. The tree only partially occludes 
the building. The occluded region contains a significant amount of information 
which we would like to use in the reconstruction process. If we are matching 
individual pixels this is simple: either a pixel is occluded and does not match 
well or it is not and matches. On theother hand, multi-pixel surfelscan contain 
both occluded and unoccluded pixels. 




I 






Mask 



Match 
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Mask 



Match 



Figure 3-7: Masks. 



More generally, either a pixel contains data which is representative of a 
surface that we intend to model, or it is an outlier. Thetree pixels in Figure 3-6 
are outliers. Other examples are transient objects (cars, people, etc), imaging 
aberrations (sensor noise, blooming), and specular reflections. To account for 
this we use a mask and allow individual pixels to be tagged as outliers and 
not count towards the match score. Deciding which are outliers and which are 
not is the hard part. One possible approach, shown on the left side of Figure 
3-7, is to tag every pixel which matches poorly as an outlier. Black indicates 
pixels tagged as outliers and poor matches. Clearly the pair (s*, s}) should not 
be considered a match. Less than 20% of the pixels are tagged yet a perfect 
match is obtained. This approach admits too many false positives, particularly 
in conjunction with shifts and illumination corrections. A more conservative 
approach is shown on the right side of Figure 3-7. A pixel is tagged as an 
outlier only if the match is equally poor throughout a small region (the eight 
nearest neighbors) around the corresponding pixel. This assumes that thetwo 
regions have already been aligned. Notice that the 1 pixel that which is truly 
an outlier is tagged as such and fewer pixels are inappropriately tagged as 
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outliers. In some cases, the number of outlier pixels is so high that the entire 
region should betaggedasan outlier. Equation 3.13 can be rewritten to account 
for both individual pixel and entire region outliers as follows: 



i/(Sj 



(3.14) 






C,P, • n, ) max 



M'&Sj e S^^C^Gfisj),^^)) / E M'tSSj e S, 



" s i eB 5 



Xs'Ss' 



E c * p i • ^ 



i(zi 



where 



Q = < 



P) e n? 




p5 e p? 


d(pj) > 


C{ — Pj 


s*- 7^ out 


ier 


E m\i^ g s,) 


1 J 3 





E!>0.5 



Hs'Ss* 



and M^lSj e Sj) returns if the pixel has been tagged as an outlier and 1 
otherwise. 

Equations 3.1, 3.2, 3.13, and 3.15 all match against a distinguished element 
p* or more generally s*. This assumes that s* is representative of the actual 
appearance of the underlying world surface. I n the presence of occlusion and 
image noise this is not always the case. Outliers in s* are not likely to match 
the corresponding data in any of the other images resulting in the mismatch 
penalty being applied |H ± | times. One way to mitigate this effect is to first 
estimates* from the image data and then perform the matching. With camera 
calibration error and variations in the viewing condition, as well as outliers in 
the image data, this is a difficult task. I nstead we successively set s* = s}. We 
can exhaustively test theses or if we are willing to occasionally miss a match 
we can simply try a few of the best. 



3.4 Discussion 



Thischapter introduces several powerful techniques for dealing with data that 
contains 1) camera calibration error; 2) significant variations in illumination; 
and, 3) difficult occlusion and image noise. As pointed out above, over-fitting 
is a concern. The next chapter applies these techniques directly to a large 
dataset. And the following chapter imposes several geometric constraints to 
prune false positives. 



Chapter 4 

From I mages to Surf els 



Chapter 3 discussed several methods to compensate for noisy data. This chap- 
ter will explore these methods in practice. We will focus on two characteristics: 
the ability to detect nearby (in both position and orientation) surfaces and once 
detected, localizetheir actual position and orientation. 

4.1 The Dataset 

A Kodak DCS 420 digital camera mounted on an instrumented platform was 
used to acquire a set of calibrated images in and around Technology Square 
(the same office complex depicted in the synthetic imagery and Figure 1-4) 
[De Couto, 1998]. Nearly 4000 images were collected from 81 node points. 
Other than avoiding inclement weather and darkness, no restrictions were 
placed on the day and time of, or weather conditions during, acquisition. The 
location of each node is shown in figure 4-1. At each node, the camera was ro- 
tated about the focal point collecting images in a hemi -spherical mosaic (Figure 
1-6). Most nodes are tiled with 47 images. The raw images are 1524 x 1012 pix- 
els and cover a field of view of 41° x 28°. Each node contains approximately 70 
million pixels. After acquisition, the images are reduced to quarter resolution 1 
and mosaiced [Coorg et al., 1998, Coorg, 1998]. Equal area projections of the 
spherical mosaic from two nodes are shown in Figure 4-2. The top node was 
acquired on an overcast day and has a distinct reddish tint. The bottom image 
was acquired on a bright clear day. Significant shadows are present in the bot- 
tom image whereas the top has fairly uniform lighting. Following mosaicing, 
the estimated camera calibrations are refined. 

After refinement, the calibration data is good, but not perfect. The pose 
estimates are within about 1° and 1 meter of the actual values 2 . As described 
in Section 3.1, these errors produce an offset between corresponding points in 
different images. A 1° pose error will displace a feature by over 8 pixels. Our 



^Si x 253 pixels. 

2 An interactive tool was used to manually identify a number of correspondences which were 
then used to bound the camera calibration error. 
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Figure 4-1: Node locations. 

calibration estimates are in an absolute coordinate frame which allows us to 
integrate images regardless of when or from what source they were collected. 
This greatly increases the quantity and quality of avail able data, but because of 
variations in illumination condition (Section 3.2) also complicates the analysis. 
Figure 4-3 and 4-4 show several images from our dataset reprqjected 3 onto a 
surfel which is coincident with an actual surface. The location, orientation, and 
size of the surfel s used are shown in Table 4.1. Surfel 1 was used to generate 
Figure 4-3 and surfel 2 Figure 4-4. If the camera calibration estimates were 
perfect and the illumination was constant, the regions in each figure should be 
identical 4 . The misalignment present in both sets is the result of error in the 
calibration estimates. Figure 4-3 is representative of the best in the dataset. A 
large number of the source images have high contrast and none of the regions 
are occluded. The third row has a distinct reddish tint. The four images in 
the center of the last row were collected under direct sunlight. And, the last 
two images were taken near dusk. Figure 4-4 is more typical of the dataset. 
It is lower in contrast and some of the regions are partially occluded by trees. 



3 U sing the estimated camera calibration. 

ignoring errors introduced during image formation (discretizing into pixels, etc) and re- 
sampling. 
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Figure 4-2: Example nodes. 



Surf el 


Position 


Normal 


Size (units) 


Size (pixels) 


X 


y 


z 


X 


y 


z 


X 


y 


X 


y 


1 


-1445 


-2600 


1200 


-1 








110 


110 


11 


11 


2 


-280 


-3078 


1500 





-1 





110 


110 


11 


11 



Table 4.1: Surfel parameters. 



62 



CHAPTER 4. FROM IMAGES TO 5URFEL5 













Figure 4-3: Reprojection ontosurfel 1 coincident with actual surface. 






Figure 4-4: Reprojection ontosurfel 2 coincident with actual surface. 
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Figure4-5: Source images for selected regions of surf el 1. 
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Figure4-6: Sou rce images for selected region s of surf el 2. 
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Figures 4-5 and 4-6 show source images with thereprojected area marked by a 
circle for several of the regions shown in Figure 4-3 and 4-4. In the worst cases 
all views of a surfel are similar to the upper left image in Figure 4-6. 

4.2 Detecting Surf els 

This section focuses on using i/(S,-) to detect nearby (both in position and ori- 
entation) surfaces. The following algorithm lists the major steps in our imple- 
mentation. 

1. Hypothesize surfel Sj in world coordinates. 

2. Select images I. 

3. Reproject regions s,-. 

4. Select nexts*. 

5. For each region s}: 

(a) Determine best shift (mj,^). 

(b) Estimate col or correction (m T) d t) m g) d g) m h) d h ). 

(c) Calculate best match. £ X {?(£$$), F(v$)) I £ 1 

6. Evaluate match set ^(S,): 

• If good enough =^> done. 

• If not =>- goto 4. 

In Chapter 2 P/s were determined by sampling along t. In essence, the im- 
age data directly determined the test points. 500 j'sand 25 a's were evaluated 
for each p*. Given the dataset described above, this would require evaluating 
~ 5 x io 12 points. Since our calibration estimates are oriented in an absolute 
coordinate system, a more efficient approach would be to choose test points in 
world coordinates. Sampling the volume of interest (from ground to 2000 units 
and containing all nodes) every 100 units in position and at 8 different orien- 
tations would test less than 2 x 10 6 points - morethan six orders of magnitude 
fewer. 

Once a surfel (such as shown in Table 4.1) is selected, we construct I. Only 
cameras which image the front side of the surfel and are not too close or too 
far are considered. We use 120 units as the minimum distance and 6000 as 
the maximum. In addition, for most of the results presented in this thesis we 
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limit |I| to 20 images. If necessary, |I| is reduced by choosing the subset which 

— _^ 

maximizes the variation in viewpoint. We use the projection ofS^Q onto the 
plane containing the surfel as a measure of the viewpoint. The subset which 
maximizes the minimum distance between the projections also maximizes the 
viewpoint variation. 

Each image in I is then reprojected on the the surfel, producing sets of 
regions Sj similar to those shown in Figures 4-3 and 4-4. To facilitate choosing 
s*, we define a region's interest as 



interest 



/ E *(^(&J),.FG5sJ))\ 



CiPj ■ a 



V~^ x&j£&_ 



V 



Vu 2 + v 2 y, 1 






Ei 



(4.1) 



where 



(u, v) e {(-l, -l), (o, -l), (l, -l), (-1, o), (l, o), (-1, l), (o, l), (l, l)} 



and then choose the largest one first. This gives priority to regions with in- 
teresting textures (i.e. ones which would produce a significant match), better 
contrast, and unforeshortened views. Only regions with a minimum interest 
(we typically use 50 as the threshold) need even be considered. In fact, at most 
we try the top fives*'s. Region 5 (top row, fifth from the left) is the most inter- 
esting in Figure 4-3 with a score of 575.5. Region 2 is the most interesting in 
Figure 4-4 with a score of 34.5 s . 

Next we consider finding the best match between s* and s*. We do this by 
evaluating: 

(4.2) 



max e[u,v 

u.v ' 



where 



e(u,v) 



x bj£bj 



xmixi^fiis*)) 






(4.3) 



Steps 5a, 5b, and 5c are actually accomplished simultaneously when finding 
the best match but for our discussion we will consider the steps separately. 

F igures 4-7 and 4-8 show -e versus u, x shift, and v, y shift, for three regions 
each from Figures 4-3 and 4-4. The left column shows the error near the cor- 
rect shift plotted as a three-dimensional surface. Thecenter and right columns 
show the projection onto the x-error plane and y-error plane respectively. Re- 
gion five from Figure 4-3 is used as s* for all of the plots in Figure 4-7. The plots 
areperiodic because of the repetitive texture (windows on the on the building). 



5 For this example we have lowered the threshold to 25 otherwise this region would not be 
considered. 
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Figure 4-7: Shift error space near correct shift (left), projection onto the x- 
error plane (center), and y-error plane (right). The periodicity is caused by 
repetitive texture (windows) on the buildings. The change of periodicity in the 
x and y directions is produced by the window spacing in the horizontal and 
vertical directions. The change in periodicity from row to row is the result of 
foreshortening. 
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Figure 4-8: Shift error space near correct shift (left), projection onto the x- 
error plane (center), and y-error plane (right). The periodicity is caused by 
repetitive texture (windows) on the buildings. The change of periodicity in the 
x and y directions is produced by the window spacing in the horizontal and 
vertical directions. The change in periodicity from row to row is the result of 
foreshortening. 
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The different periodicity in there and y directions reflects the horizontal and 
vertical spacing of the windows. The first, ninth, and last regions of Figure 4-3 
were used as s} for the first, second, and third rows of Figure 4-7 respectively. 
Region two from Figure 4-4 is used as s* for all of the plots in Figure 4-8. The 
first, fourth, and eleventh regions of Figure 4-4 were used as sj for the first, 
second, and third rows of Figure 4-8 respectively. 

The sj's used in the first and third rows of Figure 4-7 and the second row 
of Figure4-8 originate from cameras which imaged the surfel froma relatively 
oblique viewpoint. Thes* for the third line of Figure 4-7 is a very low contrast 
image acquired at dusk andthes* forthefirst lineof Figu re 4-8 contains signifi- 
cant soft occlusion from tree branches. Notethata shift is required in each case 
to obtain the best match. All of the plots are smooth with well defined minima 
making Equation 4.2 a good candidatefor optimization techniquestPressetal., 
1992]. Wehavetried a number techniques including direction set methods (eg. 
Powell's method) and conjugate gradient methods 6 and all worked well. Mul- 
tiple minima, however, area concern. We consider only the nearest minimum 
and use S i to limit the range of valid shifts. For the results presented in this 
thesis we limited u and v to ±5. In many cases, only a single minimum exists 
within this limited range. The multiple minima, particularly those in Figure 
4-7, can lead to false positives I ike the one illustrated in Figure 3-3. 

The color correction is determined using linear regression [Press et al., 
1992] and limiting the correction, sj is used as the x data and s* as the y. 
Each color channel (r, g, and b) is computed separately. Changes in both illu- 
mination and reflectance are factored into equation 3.10, however we assume 
that changes in illumination dominate the correction. Given this assumption, 
m r , m g , and m b are constrained to be positive. I mages collected under a clear 
sky tend to have a blue tint and those collected under cloudy skies have a red- 
dish tint. While changes in the spectral composition of natural lighting are 
significant they are limited. Consider the vector [m r ,m g ,m h \; If brightness is 
theonly change between s) and s*, the vector will be in thedirection [1, 1, 1]. We 
use the angle between [m r ,m g ,m h } and [1,1,1] as a measure of the variation in 
spectral composition and limit it to 10°. Further, we limit the overall change in 
brightness for any channel (m) to a maximum of 10 and a minimum of 1/10. If 
the best correction exceeds any of these limits, the match is not allowed. 

Figures 4-9 and 4-10 show the regions from Figures 4-3 and 4-4 with the 
shifts and color corrections applied 7 . Once^,^) and (m T ,d r ,m g ,d g ,m h ,d h ) have 
been determined calculating the match score is a straight forward evaluation 



6 Deriving V„,„e(M, v) is straight forward and depends only upon the image data and the 
gradient of the image data. See Appendix B for details. 

7 The correction for the 4 th region in row 3 of Figures 4-3 and 4-9 which corrects a specular 
reflection in the window exceeds the limits described above. 
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Figure 4-9: Aligned and corrected regions for surfel 1. 
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Figure 4-10: Aligned and corrected regions for surfel 2. 
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Region 


Match 


Score 


Unique 


Shift 


Slope 


1 ntercept 


X 


y 


r 


g 


b 


r 


g 


b 


1 


X 


-14.8 


29.2 


-2.5 


1.7 


1.5 


1.5 


1.6 


-29.1 


-28.5 


-47.2 


2 


X 


-15.0 


28.6 


-2.3 


1.6 


1.5 


1.5 


1.6 


-25.8 


-28.3 


-38.8 


3 


X 


-24.4 


18.3 


-2.0 


-1.9 


1.4 


1.5 


1.5 


-35.2 


-33.7 


-48.4 


4 


X 


-20.3 


22.1 


-2.0 


-1.1 


1.5 


1.5 


1.5 


-39.5 


-33.9 


-46.2 


5 


X 


-0.0 


oo 


0.0 


0.0 


1.0 


1.0 


1.0 


0.0 


0.0 


0.0 


6 


X 


-6.8 


64.0 


-0.0 


0.0 


1.0 


1.0 


0.9 


-5.3 


-0.5 


4.3 


7 


X 


-31.8 


14.9 


-3.3 


-0.7 


1.3 


1.4 


1.5 


-15.4 


-12.5 


-37.0 


8 


X 


-23.7 


19.8 


-3.2 


-0.5 


1.3 


1.4 


1.4 


-4.5 


-11.9 


-14.1 


9 


X 


-15.1 


30.0 


-1.6 


-0.6 


1.3 


1.4 


1.5 


-17.4 


-19.4 


-38.2 


10 


X 


-25.4 


18.4 


-1.6 


0.1 


1.3 


1.4 


1.4 


-14.6 


-18.8 


-23.6 


11 


X 


-34.5 


13.9 


-3.9 


1.4 


1.4 


1.4 


1.6 


-31.2 


-24.7 


-50.2 


12 


X 


-40.6 


12.0 


-4.3 


0.8 


1.5 


1.5 


1.5 


-33.1 


-32.9 


-36.7 


13 


X 


-50.7 


9.7 


-5.0 


1.4 


1.2 


1.3 


1.4 


-9.1 


-16.3 


-21.4 


14 




-152.4 


3.8 


-4.6 


-0.8 


3.7 


3.9 


3.9 


-194.0 


-153.0 


-222.8 


15 




-132.8 


4.1 


-5.2 


0.2 


4.2 


4.0 


4.1 


-245.1 


-156.9 


-226.1 


16 




-329.2 


1.7 


-1.5 


-0.0 


5.4 


5.3 


5.1 


-569.0 


-384.2 


-491.7 


17 




-1409.6 


0.9 


2.3 


0.2 


-1.2 


-1.3 


-1.4 


232.4 


195.2 


256.5 


18 




-781.6 


1.1 


7.0 


-3.4 


2.0 


1.7 


1.6 


-176.6 


-86.9 


-102.0 


19 




-162.4 


3.3 


-7.0 


2.5 


3.8 


3.9 


3.7 


-190.0 


-140.3 


-224.3 


20 




-190.0 


3.1 


-7.0 


1.1 


3.5 


3.6 


3.4 


-164.3 


-119.6 


-198.7 


21 




-132.0 


4.0 


-7.0 


0.9 


4.0 


3.9 


3.7 


-216.6 


-137.1 


-218.0 


22 




-128.3 


4.2 


-6.6 


-1.1 


1.4 


1.5 


1.7 


1.7 


-5.9 


-31.8 


23 




-904.4 


1.1 


-5.2 


-2.5 


1.1 


0.9 


0.7 


-104.0 


-55.6 


-19.9 


24 




-960.6 


1.1 


-2.9 


-1.0 


0.9 


0.9 


0.5 


-74.4 


-54.5 


-4.2 


25 




-1115.9 


1.0 


-4.0 


4.2 


0.7 


0.5 


0.3 


-38.2 


-0.8 


48.2 


26 




-924.5 


1.1 


-1.0 


-0.4 


0.8 


0.9 


0.7 


-72.9 


-57.7 


-32.2 


27 




-1052.3 


1.0 


-3.8 


-3.0 


5.4 


6.2 


0.0 


-181.7 


-128.3 


96.1 


28 


X 


-66.4 


6.8 


-0.7 


-0.6 


4.8 


5.4 


4.6 


-153.1 


-107.5 


-187.8 



Table 4.2: Match data for surfel 1. 



Region 


Match 


Score 


Unique 


Shift 


Slope 


1 ntercept 


X 


y 


r 


g 


b 


r 


g 


b 


1 




-177.2 


0.9 


-1.4 


6.0 


1.0 


1.0 


0.9 


55.5 


49.2 


73.2 


2 


X 


-0.0 


oo 


0.0 


0.0 


1.0 


1.0 


1.0 


0.0 


0.0 


0.0 


3 


X 


-21.3 


2.7 


-0.1 


1.1 


1.0 


1.2 


1.2 


15.0 


-13.5 


-24.6 


4 


X 


-17.4 


3.2 


-0.1 


0.7 


1.2 


1.2 


1.2 


-3.5 


-10.9 


-14.6 


5 




-116.3 


1.0 


0.0 


-0.1 


2.1 


2.1 


1.1 


-37.2 


-47.1 


23.5 


6 




-249.6 


0.8 


-7.1 


1.9 


0.8 


0.8 


0.4 


53.1 


31.2 


90.1 


7 




-232.5 


0.9 


-7.0 


-2.1 


2.2 


0.9 


0.1 


-86.0 


21.4 


125.8 


8 




-233.2 


0.9 


0.0 


0.0 


-1.0 


-0.9 


-1.0 


213.9 


162.7 


233.4 


9 


X 


-31.7 


1.9 


1.5 


1.5 


1.5 


1.7 


1.5 


-24.7 


-32.9 


-19.3 


10 




-116.2 


1.0 


-0.2 


-0.0 


1.7 


1.7 


1.5 


-34.8 


-16.0 


-3.2 


11 




-75.3 


1.3 


-0.3 


-1.9 


1.9 


1.9 


1.8 


-37.2 


-26.2 


-18.4 


12 


X 


-41.2 


1.8 


2.1 


1.0 


1.7 


1.7 


1.5 


-40.0 


-19.4 


-16.0 


13 


X 


-44.9 


1.8 


1.2 


2.4 


1.6 


1.7 


1.6 


-14.3 


-17.4 


-13.1 


14 


X 


-59.9 


1.5 


2.4 


2.3 


1.7 


1.7 


1.5 


-34.5 


-19.5 


-10.8 


15 


X 


-46.5 


1.8 


1.9 


0.2 


1.3 


1.5 


1.3 


-5.6 


-15.6 


-5.2 


16 


X 


-44.1 


1.8 


0.5 


0.2 


1.4 


1.6 


1.5 


-7.1 


-21.1 


-17.6 


17 




-79.8 


1.2 


-2.2 


0.0 


4.5 


4.3 


5.4 


-56.7 


-0.8 


-100.3 


18 




-69.1 


1.2 


-1.8 


0.8 


4.5 


5.3 


4.4 


-50.3 


-21.7 


-65.0 



Table 4.3: Match data for surfel 2. 
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(y+vi s i; 



^S s j)) / E 1- To help determine the significance of the 



x dj Ebj 



of E X{?Q**rj>~ 

x Oj^zOj 

match between sj and s* we use a measure called uniqueness which is related 
to the sharpness of the correlation function peak. A match's uniqueness is the 
rati oof the best match to the average match of the eight nearest neighbors and 
is defined to be 



uniqueness 





A 1 

/ u,v 


u,v Vu 2 +v' 2 E 1 



(4.4) 



jpbjfcbj 



^(^:^)^(^))/ e i 



x bj G bj 



where 



(«, u) e {(-1, -l), (0, -l), (l, -l), (-1, 0), (l, 0), (-1, l), (0, l), (l, l)}. 

A region is considered a match iftheshift and color correction arewithin limits, 
at least half of the points in s* contribute to the match and the match score is 
good enough and unique enough. For the results presented in this thesis we 
require the match score to be > -100 and the uniqueness 8 to be > 2. Tables 4.2 
and 4.3 show the match data for the regions in Figures 4-9 and 4-10. 

Finally, we evaluate the match set. We retain a surfel if the match set 
contains a minimum number of regions, the regions come from a minimum 
number of nodes, v(Sj) is greater than or equal toa minimum value and s* has 
an interest which is greater than or equal toa minimum value. For the results 
presented in thisthesis we requirethe number of regions to be > 6, thenumber 
of nodes to be > 5, z/(S.,) > -100, and interest > 50. The sets of regions shown in 
Figures 4-3 and 4-4 both produce valid matches. 



4.2.1 Results 

To test the detection characteristics of our algorithm, we evaluated z/(S,) for 
many surfels near an actual surface. A 100 x 100 x 200 unit volume, theshaded 
area in Figure 4-11, was selected so that an actual surface divides the surface 
into two cubes. 396 positions chosen at 20 units intervals were used as test 
points. Each test point was evaluated every 5° for azimuths ±45° and every 
5° for elevations ±45°. A total of 142,956 surfels were evaluated. Figure 4- 
12 shows the fraction of the nearly 13,000 surfels tested at each displacement 
that produced a valid match set. The left side of Figure 4-13 shows the frac- 
tion of detections versus displacement and the angle between the actual and 



8 While growing surfaces (discussed in the next chapter) we allow uniqueness as low as 1.5. 
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Figure 4-11: Test volume (small Figure 4- 12: Average detect ion rate, 
shaded rectangle). 
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Figure 4-13: Detection rates. 
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Figure 4-14: Empty test volume 
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Figure 4-15: Scale effects for an 
un occluded surface. 
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Figure4-16: Scale effects for an occluded surface without (left) and with (right) 
masks. 
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tested orientation. The right side shows the fraction of detections versus az- 
imuth and elevation 9 . We also investigated theeffect of surfel size on detection 
rate. As a control, we measured the detections in the empty volume shown 
as a shaded region in Figure 4-14. An actual surface was selected and tessel- 
lated with surfels. The fraction of false negatives (1- the detection rate) for 
an unoccluded surface versus the fraction of true negatives (l- the detection 
rate) for the empty volume is shown in Figure 4-15. Surfel sizes from 5 units 
(lower left corner) to 50 (upper right corner) in increments of 5 units are plotted 
as diamonds. The "bump" in Figure 4-15 is the result of noise and the discrete 
nature of surfels. Similar curves for an occluded (significant tree coverage) sur- 
face are shown in Figure 4-16. The left and right plots in Figure 4-16 a re for 
the same surface. The data shown in the right plot utilizes the mask function 
described in Section 3.3. A surfel size near 10 units gives the best all around 
performance. The best surfel size is a function of the data and is related to 
the size in world coordinates of the significant image features- in our casethe 
windows. Potentially scale-space techniques such as [Li and Chen, 1995] can 
be used to determine the best surfel size. 



4.3 Localizing Position and Orientation 

As shown in the previous section, our method is capable of detecting surfaces 
a significant distance (both in position and orientation) from S,. The algorithm 
described so far produces a set of corresponding textures, %. The underlying 
surface which gave rise to the textures may have a position and orientation dif- 
ferent from that of S.,. The i nformation contai ned i n s,- can be used to esti mate 
the position and orientation of the underlying surface. 

Once a surface has been detected, we localize its position and orientation 
using the foil owing algorithm: 

1. Until convergence: 

(a) Update surfel position P,. 

(b) Update surfel orientation n,-. 

(c) Reevaluate //(Sj). 

Steps la, lb, and lc are actually interdependent and ideally should be calcu- 
lated simultaneously. For ease of implementation and to speed convergence 10 



9 AII surfels tested at a particular azimuth and elevation regardless of displacement are 
included in the fraction. 

10 Theoretically, convergence is 0(n 3 ) where n is the number of parameters optimized, thus 
two half-sized problems converge four times faster than the full-sized one. We have observed 
this speed up in practice. I n addition, the symbolic gradient expression needed by conjugate 
gradient methods is much simpler for the two half-sized problems. 
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we have separated them. As will be seen in this section, we obtain good results 
performing the steps individually. 

We use a slightly modified version of Equation 3.4, which considers only 
matching regions, to update P.,: 



argmin ]T Qp} x ((P,- - Q) x Qpj) • (P,- - Q). (4.5) 

When Pj is updated, u { and ^ must also be updated, s* changes with the new 
Pj but p* does not, thus (w», ij) must be corrected for the difference between the 
old and news*. 

A match's dependence upon S,- andn, can be made more explicit by rewriting 
Equation 3.15 as follows: 



where 



^fCiPj-nj-Jmaxefo) 
u(Sj) = ^ =^- (4.6) 

ieQ 

Y, xwrQPi,?) + M,.f(t(^,i*))) 



v^ 



r(u ; -) - ^^ =^ . (4.7) 



y 



x Oj £= kj 



To update n,- we evaluate: 

argmax V] e(n.,). (4.8) 

Both direction set methods and conjugate gradient methods 11 work well to 
solve Equation 4.8. 

Once the position and orientation of thesurfel have been updated, wereper- 
formSteps2, 3, 5, and 6 of the detection algorithm. Step 2 is only partially per- 
formed. Images with currently matching regions are retained. Images which 
no longer view the front side of thesurfel are eliminated and additional ones 
that do are added. Steps la, lb, and lc of the localization algorithm are re- 
peated whilethe match score is improving until the position update is < 1 unit 
and the orientation update is < 1° or for 3 iterations if the match score is not 
improving. 



"•While not as straight forward as in the last section, deriving V nj e(n.,-) leads to an easily 
evaluated expression that depends only upon known quantities such as the image data, the 
gradient of the image data, and the position and orientation of surfel. See Appendix B for 
details. 
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Figure4-17: Surfel localization. 
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Figure 4-18: Localization examples. 
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Figure 4-19: Localization summary. 
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4.3.1 Results 

Figure 4-17 depicts the localization process. I mages which view thefront side of 
thesurfel are shown in grey. Those that contribute to the match set are shown 
in black. Lines of sight through the center of the matched regions and building 
outlines are also shown in grey. Figure 4-17a shows the initial detection for 
a surfel displaced from the actual surface by 100 units in positions and 3CP 
in orientation. Figure 4-17b shows the surfel after updating its position and 
orientation (Steps la and lb of the localization algorithm). F igu re 4-17c shows 
a new match set using the updated position and orientation (Steplc). Figure4- 
17d shows the final match set at convergence. The same set of test points used 
to test detection (Figures 4-12 and 4-13) was also used to evaluate localization. 
Figure 4-18 shows two test points evaluated at 361 different orientations. The 
test point for the plot on the left is 80 units in front of the actual surface and 
the one for the right is 100 units behind. The asymmetry in both of these plots 
is caused by the asymmetric distribution of cameras shown in Figure 4-24. A 
diamond is plotted at the initial orientation for each detection and a pi us marks 
the final estimated orientation. Figure 4-19 shows the aggregate results for 
the complete set of test points. The plot on the left shows final displacement 
versus the initial displacement and the angle between the actual and tested 
orientation. The plot on the right shows the angle between the actual and final 
estimated orientation versus the initial displacement and the angle between 
the actual and tested orientation. For displacements upto 100 units in position 
and 30° in orientation, nearly all of the test points converged to the correct 
values. 



4.4 Discussion 

So far we have focused on detecting and localizing nearby surfels and have not 
addressed bogus matches. As pointed in Chapter 3, compensating for noisy 
data admits false positives. Figure 4-20 shows the raw surfels 12 detected and 
localized near one of the buildings imaged in thedataset. There are a signif- 
icant number of false positives. Notice the parallel layers of surfels near the 
upper right corner of the building. These are false positives similar to the one 
shown in Figure 3-3. Many of the surfels are near the surfaces of the building. 
Figure 4-21 shows the distribution of distances to the nearest model surface. 
Figures 4-22 and 4-23 show the raw results for the full reconstruction and Fig- 
ure 4-24 shows the volume used for the full reconstruction. The next chapter 
addresses ways to el i mi nate fal se positives. 

The reconstructions shown in this thesis were performed on quarter reso- 



12 Only front facing surfels are rendered. 
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Figure 4-20: Raw surfels (partial reconstruction). 
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Figure 4-21: Distance to nearest model surface (partial reconstruction). 
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Figure4-22: Raw surf els (full reconstruction). 
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Figure4-23: Distance to nearest model surface (full reconstruction). 
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Figure4-24: Reconstruction volume (shaded area) used for full reconstruction. 
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Figure 4-25: Raw surfels for full reconstruction using half (top) and eighth 
(bottom) resolution images. 
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lution images. To test the effects of image resolution, we reran the raw detec- 
tion and localization algorithms on half and eighth resolution images. Several 
parameters (maximum shifts, interest, and uniqueness) are resolution depen- 
dent and these were adjusted appropriately. The results are shown in Figure 
4-25. The half resolution reconstruction is very sparse. We suspect the primary 
cause is image noise. It turns out that half theoutput resolution actually corre- 
sponds to the full resolution of the physical sensor [DeCouto, 1998]. Increased 
image noise degrades the match score, whose tolerance was not adjusted, and 
interferes with the gradient estimates used to optimize the shifts. Assuming 
the noise is Gaussian in nature, reducing the resolution also reduces the noise. 
However, there is a limit to how much the resolution can be reduced; sufficient 
information must be retained to perform the reconstruction. This implies that 
there is a range of image resolutions (with favorable noise characteristics and 
sufficient information content) which are essentially equivalent. The quarter 
and eighth resolution reconstructions seem to support this view. 



Chapter 5 

From Surfelsto Surfaces 



Theshifts, illumination corrections, and masks introduced in Chapter 3tocom- 
pensate for noisy data also increase the occurrence of false positives. The re- 
sults presented in Chapter 4 are purely local and make no attempt to reject 
these false positives. This chapter explores several geometric constraints which 
together eliminate nearly all false positives. 

5.1 Camera Updates 

As poi nted out i n Section 3.1 unconstrained shifts introduce false positives such 
as the one shown in Figure 3-3. Figures 5-1 and 5-2 shows shifts {(ii,v)} plot- 
ted as a vector field in image-space along with the corresponding image. The 
position of each vector corresponds to the location of the center of a matching 
region. Note that theshifts appear to be random. The actual image displace- 
ments caused by camera calibration error should be consistent with a transla- 
tion of the camera center and a rotation about it or formally 



where 

and Pj is the location of S 
algorithm: 

1. For each camera P. 

(a) Build {(ui,Vi)^}. 

(b) Calculatea first-order camera calibration update,?'. 

(c) Calculate thefinal camera calibration update,!*". 

(d) Remove matching regions with shifts that are not consistent with the 
final update, P". 
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To enforce this constraint we 
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the foil 


owing 
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Figure 5-1: Shifts {(ui,Vi)} plotted as a vector field and image data for node 27 
image 13. 
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Figure 5-2: Shifts {{u h Vi)} plotted as a vector field and i mage data for for node 
28 image 17. 
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Figure 5-3: Effects of camera calibration updates. 



2. Prune surf els which no longer meet the match criteria. 

The left side of Figure 5-3 shows the effect of translating a camera in thef 
direction. Ax and TV are related as 



Ax 



f 



CP-z 



(5.2) 



and similarly for Ay and T s 



Ay 



f 



CP-z 



-T- 

±y. 



(5.3) 



The right side of Figure 5-3 shows rotation about they axis. Differentiating 

' X — Xq\ 



9 = arctan 



yields the first order approximation 



./ 



Ax=- f2+ (y^ 2 se. 



(5.4) 



Similarly for </>, rotation about there axis, 



Av = /Lt^%. 



(5.5) 
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Finally, for 7, rotation about the optical axis, 

Ax = -(y-y )5-f (5.6) 

Ay = (x-x )5-f. (5.7) 

Equations 5.2-5.7 combi ne to form the following pair of coupled linear equa- 
tions: 

u = J - {X ~ Xo) 89 + =j— T* + (y - y )6>y (5.8) 

(x - x )8 1 . (5.9) 

J CP-z 

Equations 5.8 and 5.9 aresolved in a least squares fashion [Watkins, 1991]. 
P is updated using 89, 84>, 5j, T s , and T/ to producethe initial solution in Step 
lb, P'. Thefinal update is 

argmin £ \\(u, v) - T(P jf V") + T(P,, P)|| 2 , (5.10) 

where 



f 2 


+ (x- 


-x y 
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-59 + 
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+ {y- 


-yo? 


-8(1) + 
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T- 



Ji, 

(ui,vi)€i 



Q= {(ii,ii) \e> \\(u,v) -T(P,,P') + T(P,-,r)l| 2 } . 

Typically e is 1 pixel. Solutions which have enough data (at least 4 points) and 
fit well are retained. P" is used to prune inconsistent matches from s,-, 

B5\i i ifc<||(ti,t;)-T(P i ,r ff ) + T(P i ,r)|| 2 . 

After removing inconsistent matches, if % no longer meets the criteria de- 
scribed in Section 4.2, it is discarded. Figure 5-6 shows the average rms shift 
before and after performing the calibration update. The number of consistent 
data points (shifts) are plotted along there axis and the square root of the mean 
squared shift is plot along the y axis. Note that the necessary shifts are sig- 
nificantly reduces. Ideally, after updating the camera calibration estimates, no 
shifts would be required. The images used as input contain a small amount of 
nonlinear distortion which cannot be modeled by the pin-hole camera model. 
We estimate the rms value of this distortion to be about 1 pixel. 

Figure 5-4 shows the consistent surfels remaining after applying this algo- 
rithm tothe raw reconstruction shown in Figure4-22 and Figure 5-5 showsthe 
distribution of distances tothe nearest model surface. A number of theconsis- 
tent surf els come from objects which are not in the reference model. The cluster 
of surfels between the building outlines near the top center of Figure 5-4 is one 
example. These surf els come from a nearby building. Asa test, a small portion 
of thefull reconstruction was rerun with the updated camera calibration. The 
results are shown in Figure 5-7. For comparison, the consistent surfels from 
the same area are shown in Figure 5-8. 



^he internal parameters are held fixed and for stability reasons TV is not updated. 
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Figure 5-4: Consistent surf els. 
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Figure 5-5: Distribution of error for consistent surfels. 




i 



\K' 



Initial calibration 
Final calibration 




500 1000 

Data points 



1500 



Figure 5-6: Comparison of initial and final calibration estimates. 
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Figure 5-7: Close-up of reconstruction using updated camera calibrations. 




Figure 5-8: Close-up of same area showing only consistent surfels from original 
reconstruction (estimated camera calibrations). 
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5.2 One Pixel One Surfel 
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Figure 5-9: A region from one image 
which contributes to multiple surfel s. 



Figure 5-10: Determining if S b 
neighbor of S Q . 



is a 



Figure 5-9 shows a region from one image which contributes to the match 
set of multiple surfel s. Each pixel in each image should contri bute to at most 
one surfel. Deciding which surfel is the hard part. Detection and localization 
as described in Chapter 4 do not enforce this constraint and as a result even 
after enforcing a consistent calibration update there are many image regions 
which contributetomultiplesurfels. We eliminate them in a soft manner using 
thefollowing algorithm: 

1. Score each surfel. 

2. For each surfel S a . 

(a) For each region s*gs„. 

i. For each surfel S 6 with a score higher than S ai if s\ e s 6 . 
A. De-weights 4 a . 

(b) If the match score is no longer sufficient, prune S a . 

To score a surfel we use a combination of the number of cameras which 
contribute to a surfel and the number of neighbors that it has. A surfel S a is 
considered a neighbor of S b if 1) the distance from P b to the plane containing 
S a (the normal distance) is no more than 4; 2) the distance from the projection 
of S b onto the plane containing S a and P a (the tangential distance) is no more 
than d t ; and, 3) the angle between n a and n 6 is no more than f3. This notion of 
neighbors is essentially a smoothness constraint and will be used in the next 
section to group surf els. Typically, we use 4 = 15, d t = 300, and/5 = arccos(0.9). 
The de-weighting is done in a continuous manner. S a is divided into sample 
points (we use 10 x 10) each with an initial weight of 1.0. Each sample point 
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Figure 5-11: Surfels after pruning multiple contributions. 
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Figure 5-12: Distribution of error for pruned surfels. 

is tested against the cone formed by S b and Q. Those inside the cone have 
their weights updated. The new weight is based on the normal distance and 
the angle between n a and n 6 . If the angle is greater then /3, the new weight 
is 0.0. Otherwise, the new weight is based on the normal distance; 1.0 if it is 
less than d z , 0.0 if it is greater than 34, and varies linearly in between. The 
lesser of the currently assigned and new weight is retained. After all $/s have 
been tested against S a , the average sample point weight is used as the new 
weight for sj,. When each s^ has been tested the sum of the weights is used to 
determine if S a is retained. Figure 5-11 shows the reconstruction after pruning 
multiple contributions and Figure 5-12 shows the distribution of distances to 
the nearest model surface. 



5.3 Grouping Surfels 



The buildings we are trying to model are much larger than an individual surfel. 
Therefore, a large number of surfels should be reconstructed for each actual 
surface. Using the notion of neighbors described in the last section, we group 
the reconstructed surfels as follows: 

1. For each surfel S a . 
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(a) For each surfel S 6 already assigned a group. 

i. If S a and S fe are neighbors. 

A. If S a has not already been assigned to a group, then assign S a 
to the group contai ni ng S 6 . 

B. Otherwise merge the groups containing S a and S b . 

(b) If S a has not been assigned to a group, then create a new group and 
assign S a to it. 

In practice we retain only groups which have at least a minimum number 
of (typically four) surfels. All of the surfels in a group should come from the 
same surface. This notion of grouping places no restriction on the underlying 
surface other than smoothness (eg. it may contain compound curves). Figure 
5-13 shows the reconstruction after grouping and removing groups with fewer 
than four surfels. Nearly all of the surfaces in the reference model have at 
least one corresponding group. Figure 5-14 shows the distribution of distances 
tothe nearest model surface. 

5.4 G rowi ng Su rf aces 

M any of the groups shown i n figure 5-13 do not completely cover the underlyi ng 
surface. There are several reasons why surfels corresponding to actual surfaces 
might not produce a valid match set. The main one is soft occlusion described 
in Section 3.3. Another is local maxima encountered while finding the best 
shifts and updating the surfel's normal. I n the reconstruction process so far, 
the mask technique described in Section 3.3 has not been utilized. In this 
section we make use of it. We also use estimated the shifts and illumination 
corrections to help place us closer to the correct maxima. We use the foil owing 
algorithm to grow surfaces: 

1. For each group. 

(a) Create an empty list of hypothesized surfels. 

(b) For each hypothesized surfel. 

i. Test using the detection and localization algorithms, 
ii. If a match. 

A. Add tothe group. 

B. Test against each surfel in each of the other groups. 
If a neighbor, merge the two groups. 

(c) Use the next surfel in the group S a to generate new hypotheses 
and goto Step lb. 



5.4. GROWING SURFACES 



97 





Figure 5-13: Surfels after grouping. 
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Figure 5-14: Distribution of error for grouped surfels. 

The hypotheses in Step lc are generated from S a by considering the eight 
nearest neighbors in the plane containing S„ as shown in Figure 5-15. S„ is 
shown in grey and the hypotheses are shown in white. Each of the other recon- 
structed surfels in the group are orthographical ly projected onto the hypothe- 
sized surfels and hypotheses which are more than half covered are discarded. 
The shifts and illumination corrections associated with S a are used as initial 
values for each hypothesis in Step l(b)i. In addition we lower the minimum 
interest to 25 and the minimum uniqueness to 1.2. Figure 5-16 shows the 
reconstruction after growing. After growing, the coverage of each surface is 
nearly complete. Figure 5-17 shows the distribution of distances to the nearest 
model surface. Figure 5-18 shows grown surfels after removing multiple con- 
tributions and grouping as described in the previous two sections and Figure 
5-19 shows the distribution of distances to the nearest model surface. 



5.5 Extracting Models and Textures 

So far, the only assumption we have made about the structure of the world is 
that locally it can be approximated by a plane. All of the buildings imaged in 
our dataset are composed of planar faces, therefore we simply fit planes to the 
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Figure 5-15: Reconstructed surfel (grey) and hypotheses (white). 

groups identified in the previous section. In this case, a face is equivalent to a 
large surfel. The position and orientation of S g are determined by 



E "P; / E 



w 



(5.11) 



and 



IX, 



E 



wn 




(5.12) 



where w is the score calculated in Step 1 of the algorithm to remove multiple 
contributions and S g is the set of all surfels in group g. The extent of S g is 
determined by orthographical I y projecting^ ontoS s and finding the bounding 
box. Figure 5-20 shows the reconstructed S/s. A total of 15 surfaces were 
recovered. Figure 5-21 shows the distribution of distances to the nearest model 
surface. As noted previously, many surfels come from structures not in the 
reference model. Three of the reconstructed surfaces fall into this category, 
hence Figure 5-21 has a maximum of 80. 

I mage data from contributors to s a can easily be related using the illumina- 
tion corrections calculated during detection and localization. The relationship 
between s\ and s a where s a and s b are members of the same group is more 
difficult when i is not a contributer of s a . To resolve these relationships we 
construct a table of the illumination corrections for J, the set of images con- 
tributing toat least one surfel in S, using the foil owing algorithm: 

1. For each group. 

(a) Find the surfel with the highest score, S m . 

i. Makes^ the root of the tree, 
ii. Add the illumination correction for s^ to the table where i ^ *. 
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Figure 5-16: Surfels after growing. 
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Figure 5-17: Distribution of error for grown surfels. 

(b) Until no new entries to the table are added, 
i. For each S.,. 

A. If 3 a, b | s" is in the table and s^ is not, 

Then calculate the illumination correction between^ and 
and enter it in the table. 



3 



We assume that each surfel in a group has at least one contributing image in 
common with at least one other surfel in the group. When building the table 
we consider only the first correction encountered for each contributing image. 
Oncethe table is built we correct each image in J and reproject it ontoS 5 . The 
texture associated with S g is simply the average. Figure 5-22 shows two views 
of the reconstructed textures. Notice that the rows of window in adjacent faces 
are properly aligned. This occurs even though no constraints between faces are 
imposed. 



5.6 Discussion 

This chapter uses several simple geometric constraints to remove virtually all 
false positives from the purely local reconstruction described in Chapter 4. Af- 
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Figure 5-18: Surfels after regrouping. 
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Figure 5-19: Distribution of error for regrouped surfels. 
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ter imposing consistent calibration updates, removing multiple contributions 
and grouping, the remaining surfels are excellent seeds for growing surfaces. 
Of the 16 surfaces in the reference model, 12 were recovered. All of the re- 
maining surfaces are severely occluded by trees. Nearly all of the images are 
similar to the upper left-hand image of Figures 4-4 and 4-6. In spite of this 
several surfels were recovered on two of the surfaces, however they did not 
survive the grouping process. Figure 5-23 shows the results of growing these 
surfels. The top shows the raw surfels after growing and the bottom shows re- 
constructed textures. I n addition to being severly occluded by trees, the other 
two surfaces have very little texture and one of them suffers from a lack of 
data. Three surfaces from adjacent buildings not contained in the model were 
also recovered. Thefacenear the top center of the upper image in Figure 5-22 
is from the Parson's lab. The surfaces on the left of the upper and the right of 
the lower image is from Draper lab. 

A summary of the computation required is presented in Table 5.1. The ma- 
jor steps of our algorithm and the section which describes them are listed along 
with the elapsed and cpu runtimes, number of surfels output, and complexity. 
A total of 64.5 hours of cpu time was needed to produce the textured model 
shown in Figure 5-22. This is less than 1 minute of cpu time per image and 
less than 0.2 seconds of cpu time per test point. The "Detection & Local iza- 
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Figure 5-20: Raw model surfaces. 
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Figure 5-21: Distribution of error for model surfaces. 



Step 


Section 


Time 


Surf els 


Complexity 


Elapsed 


cpu 


Detection & 
Localization 


4.2 
4.3 


42hrst 


60hrs 


54,212 


0(V xNx S) 


Camera 
U pdates 


5.1 


18.5min 


17min 


2977 


0(N x #?) 


1 Pixel, 
1 Surfel 


5.2 


3.5min 


2min 


1636 


0{N x # 2 ) 


Group 


5.3 


4min 


2.5min 


1272 


0(# 2 ) 


Grow 


5.4 


6hrst 


4hrs 


3007 


0(A x N) 


Model 


5.5 


11.25min 


9.5min 


15 


0(#) 


Texture 


5.5 


24.5min 


15.5min 


15 


0(A x N) 



Table 5.1: Run times and orders of growth. 
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Figure 5-22: Textured model surfaces. 
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Figure 5-23: Textured model surfaces with two additional surfaces. 
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tion" and "Grow" steps were performed in parallel on a cluster of 32 400Mhz 
dual processor Pentium II machines running I inux. The remaining steps were 
performed on a single 200Mhz uniprocessor Pentium Pro machine also run- 
ning I inux. The two times marked with daggers reflect the speed up of parallel 
processing. All others (including cpu time for "Detection & Localization" and 
"Grow") are total times. The complexities shown are upper bounds. V is the 
reconstruction volume, N is the total number of images, S is the surfel size 
(area in world coordinates), # is the total number of recovered surfels, and A 
is the total surface area (in world coordinates) of the scene to be reconstructed. 



Detection & Localization 

As noted above, the reported elapsed time reflects parallel execution. Thetotal 
elapsed time is 1,260 hrs (42 hrs x 30 nodes). The cpu time required is 1/20 
of the total elapsed time. The major cause of this is I/O. The reconstruction is 
divided into 6000 chunks. Each of these chunks is reconstructed independently 
and requires loading a large image set from a file server across a network. The 
work required to test a single surfel is linearly proportional to surfel area in 
world coordinates (S) and the number of images which view the surfel. The 
total number of surfels tested is linearly proportional to the reconstruction vol- 
ume. N is an upper bound on the number of images which view a surfel, but 
it is not a tight bound. For example, only a fraction of the dataset can image 
a particular surfel. Therefore, if the dataset is acquired with roughly a fixed 
node density and images only contribute to surfels which are not too far away 
(Section 4.2), we would expect the number of images which view a surfel to 
eventually become independent of TV. In essence, once the dataset has reached 
a certain size, new images extend the reconstruction volume but do not affect 
the local reconstructions at the core of the volume. For large data sets the 
expected order of growth is 0(V x S) 



Camera Updates 

The work required to update a single camera is proportional to the number of 
surfels it contributes to squared and a total of N cameras must be updated. # 
is an upper bound on the number of surfels which a camera can view, but it is 
not a tight one. Similar to the discussion in "Detection & Localization", once 
the reconstruction reaches a certain size, the number of surfels imaged by a 
given camera should remain roughly constant. Spatial hashing can be used to 
help exploit this property. Thus, for large reconstructions the expected order of 
growth \sO(N) 
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1 Pixel, 1 Surfel 

The work requi red to i impose the "1 Pixel, 1 Surfel " constrai nt for a si ngle i mage 
contributing to a given surfel is linearly proportional tothenumber of surfelsin 
a cone formed by the camera's center of projection and the surfel and extending 
out to the maximum allowed distance from the camera. This constraint must 
be i mposed for each image contributing to each of the # surfel s. If images can 
contribute only to surfels which are not too close (Section 4.2), then once the 
reconstruction reaches a certain size, the number of surfels in the cone should 
remain roughly constant. Therefore, for large reconstructions with large data 
sets the expected order of growth is 0(#). 

Group 

The work required to group a surfel is linearly proportional to the number of 
surfelsin the vicinity and all # surfels must be grouped. The number of surfels 
in a given area is bounded and by using spatial hashing the expected order of 
growth becomes 0(#). 

Grow 

The comments in "Detection & Localization" about elapsed time apply here, 
except that only 15 surfaces were grown, therefore only 15 nodes were used. 
Assuming that the growing algorithm effectively stops at surface boundaries, 
there are a total A/S locations on all faces to be tested. For the faces shown in 
Figure 5-22 this is a good assumption. For large data sets, each location tested 
requires O(S) work, giving an expected order of growth of 0(A). 

Model 

The work required to fit a (plane) model is linearly proportional tothenumber 
of surfels in the group. Since every surfel belongs to a group, the given order of 
growth 0(#) is a tight bound. 

Texture 

The work required toextract the texture associated with a face is linearly pro- 
portional to the area of the face and the number of images which view it. N 
is an upper bound on the number of images which can view a face and si nee a 
face can be of arbitrary size it is not clear that we can do better. Therefore, the 
order of growth is 0(A x TV). 
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Chapter 6 
Conclusions 



This thesis presents a novel method for automatically recovering dense surfels 
using large sets (1000's) of calibrated images taken from arbitrary positions 
within the scene. Physical instruments, such as Global Positioning System 
(GPS), inertial sensors, and inclinometers, are used to estimate the position 
and orientation of each image. Long baseline images improve the accuracy; 
short baselines and the large number of images simplify the correspondence 
problem. The initial stage of the algorithm is completely local enabling par- 
allelization and scales linearly with the number of images. Subsequent stages 
areglobal in nature, exploit geometric constraints, and scalequadratically with 
the complexity of the underlying scene. 
We descri be techniques for: 

• Detecting and localizing surfels. 

• Refining camera calibration estimates and rejecting false positive surfels. 

• Grouping surfels into surfaces. 

• Growing surfaces along a two-dimensional manifold. 

• Producing high quality, textured three-dimensional models from surfaces. 
Some of our approach's most important characteristics are: 

• It is fully automatic. 

• It uses and refines noisy calibration estimates. 

• It compensates for large variations in illumination. 

• It matches image data directly in three-dimensional space. 

• It tolerates significant soft occlusion (eg. tree branches). 

• It associates, at a fundamental level, an estimated normal (eliminating 
thefrontal-planar assumption) and texture with each surfel. 

Ill 
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Our algorithms also exploit several geometric constraints inherent in three- 
dimensional environments and scale well to large sets of images. We believe 
that these characteristics will be important for systems which automatically 
recover large-scale high-quality three-dimensional models. A set of about 4000 
calibrated images was used to test our algorithms. The results presented in 
this thesis demonstrate that they can be used for three-dimensional recon- 
struction. To our knowledge, the City Scanning project (eg. [Coorg, 1998] 
and the work presented in this thesis) is the first to produce high-quality tex- 
tured models from such large image sets. The image sets used in this thesis 
are nearly two orders of magnitude larger than the largest sets used by other 
approaches. The approach presented in this thesis, recovering dense surfels 
by matching raw image data directly in three-dimensional space, is unique 
amoung the City Scanning approaches. 

6.1 Future Work 

The major limitation of the work presented in this thesis is the difficulty of 
performing reconstruction "through the trees". A large part of the difficulty 
is caused by matching against s*. If s* is not representative of the underlying 
texture, then the match set s., will most likely be insufficient resulting in a 
false negative. If we knew the underlying texture and used it ass*, the match 
set should improve significantly. Similarly, estimating the underlying texture, 
rather than simplying selecting one of the views, and using it ass* should also 
improve the match. Potentially this would allow us recover the missing faces 
in Figure 5-22. 

Several other areas could also benefit from further exploration: 

• Develop alternate techniques for exploring the volume of interest. 

• Explore improvements/optimizations to the basic algorithms. 

• Develop a more sophisticated model of illumination and reflectance to im- 
prove the correction. 

• Explore fitting a broader class of models to surfels. 

• Develop techniques to enhance the recovered texture. 

6.1.1 Exploring the Volume of Interest 

As noted in Section 4.2, we exhaustively test all of the nearly 2 x 10 s points 
in the volume of interest. About 2% of the tested points produce surfels and 
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Figure 6-1: Free space. 

about 6% of those are consistent with a single set of camera calibration up- 
dates. One way of reducing the number of points that need to be tested is to 
use techniques such as binocular stereo to generate test points that are likely 
to produce surfels. Because of error in the estimated three-dimensional posi- 
tion, a small volume around the point should be tested. If a surfel is detected 
and localized, additional test points can be generate nearby in a manner sim- 
ilar to that used in Section 5.4 to grow surfaces. Another possible approach 
is to keep track of empty space. For example, if surfel S in Figure 6-1 corre- 
sponds to an actual surface, then the shaded area must be empty space and 
need not be tested. Conservatively estimating regions of empty space during 
the reconstruction process could significantly reduce the number of points ac- 
tually tested. Argus, our image acquisition platform is pushed for node to node 
while obtaining a dataset and the path it follows also identifies empty space 
that limit test points. I n a similar fashion, it may prove useful to keep track 
of the general location of soft occluders. For example, views which do not have 
soft occlusion could be selected preferentially or weighted more heavily. 

6.1.2 Improving the Basic Algorithm 

Several components of the basic algorithm could be improved. Several param- 
eters such as surfel size and thresholds for interest and uniqueness are se- 
lected and utilized on a global basis. These values are dependent upon the 
data and their values should reflect this dependence. For example, the overall 
image brightness and contrast havea large impact on interest and uniqueness 
scores. For various reasons (geometry, sun position, etc.), some building faces 
are frequently brightly illuminated and others are nearly always in shadow. 
Using the same threshold for both tends to produce a large number of false 
positives in theformer case and a large number of false negatives in the latter. 
Thresholds based upon the imaging conditions should produce better results. 
The quality of matches could also be improved by performing comparisons at 
the appropriate resolution. The effective resolution of reprqjected regions vary 
with viewing condition (eg. distance and foreshortening). Source images could 
be loaded into a pyramid of resolutions and prior to matching, a level selected 
such that all the reprqjected textures are of similar resolution. This technique 
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should significantly improve the algorithms ability to match distant or highly 
forshorted views. As noted in Section 4.3, updating a surfels position and 
orientation should ideally be done simultaneously. Formulations which com- 
bine these optimizations may localize surfels more accurately. And finally, as 
noted in Section 3.3, our masking algorithm is conservative. More sophisti- 
cated methods of identifying outlier pixels would enable the masks to be more 
widely used and should result in fewer false negatives. 



6.1.3 Illumination and Reflectance Models 

The correction proposed in Section 3.2 does a good job of compensating for 
changes in viewing condition (illumination, view point, reflectance, etc.), but 
theextra degrees of freedom also admit false positives such as shown in Figure 
3-5. Oneway to reduce the number of false positives is to constrain the correc- 
tions by measuring the viewing conditions. For example, the total irradiance 
arriving at a node and an estimate of cloud cover might be sufficient for a rough 
estimate of the correction. The corrections used to extract the textures shown 
in Figure 5-22 are consistent within a face, but not necessarily between faces. 
A direct estimate might make it possible to locate all corrections in a global 
correction space. 

As demonstrated in Figures 4-3 and 4-9 in some cases our method is ca- 
pable of compensating for specular reflection. These corrections are generally 
rejected in an attempt to limit false positives. One way to improve match- 
ing for surfaces with substantial specular characteristics is to decompose the 
reflectance into specular and diffuse components and use just the diffuse com- 
ponent for matching. Nayar et al. [Nayar et al., 1997] show that the specular 
component of an objects reflectance tends to be linearly polarized while the 
diffuse component is not. Using this property and a polarization image 1 , Na- 
yar et al. present an algorithm for separating the components. Wolff [Wolff 
and Andreou, 1995, Wolff, 1997] demonstrates a practical device for acquiring 
polarization images. 

Finally, given a number of images of thesame world surface taken under dif- 
ferent viewing conditions, it should be possible to estimate its bidirectional re- 
flectance distribution function (BRDF)tHorn, 1986]. Tagareand DeFigueredo 
[Tagareand deFigueiredo, 1993], Oren and Nayar [Oren and Nayar, 1995], and 
Dana et al. [Dana etal., 1996] present a techniques for fitting BRDF's to image 
data. 



X A set of color images of the same scene acquired from the same location each imaged with 
a polarization filter at a different orientation 
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6.1.4 More Descriptive Models 

The models extracted in Section 5.5 are composed of simple rectangles. The 
buildings in our dataset are well modeled by planar faces, however many build- 
ings have components which arenot. Fitting general modelstodata isan open 
problem, however expanding the set of primitives to include cones, cylinders, 
and spheres as well as planes would greatly increase the descriptive power of 
our models. Further, the faces recovered should be part of a closed solid. This 
constraint could easily be enforced by hypothesizing missing faces. Imposing 
this constraint would result in recovering two of the four un recovered faces. As 
part of the City Scanning project, Cutler [1999] has investigated aggregating 
reconstructions from multiple sources and imposing solid constraints. 

6.1.5 Better Textures 

As described in Section 5.5 the textures displayed in Figure 5-22 aregenerated 
by averaging the illumination corrected image data. The results are good but 
can be improved. Iteratively estimating the average texture should improve 
the results. For example, thecurrent average texture is used to reesti mate the 
illumination corrections for each image and weights for each pixel of each im- 
age. The new illumination corrections and weights are then used to reesti mate 
the average texture until it converges. Noattempt has been madetoalign the 
data from individual images. Warping techniques similar to those described by 
Szeliski [Szeliski, 1996, Szeliski and Shum, 1997] should improve the quality 
of the extracted texture. 
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Appendix A 
Camera Model 
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FigureA-l: Camera model. 



The pi n-hole camera model uti I ized throughout this thesis uses the perspec- 
tive projection model of image formation. A4x3 matrix maps points expressed 
in homogeneous coordinates from V 3 to V 2 . The left side of Figure A-l shows 
a camera centered Cartesian coordinate system. The camera's optical axis is 
coincident with the z axis and its center of projection (C) is at the origin. The 
image plane is parallel toXhexy plane and located a distance /from the origin. 
Even though the image plane is not required to be parallel tothexy plane, most 
camera models do not explicitly consider this possibility. I nstead the effects of 
image plane pitch 6> x and tilt 9 y are lumped in with lens distortion under the 
heading of thin prism distortion. We assume that 6> x = 9 y = for consistency 
with traditional pi n-hole camera models. The point where the image plane and 
the optical axis intersect is known as the principal point. Under perspective 
projection a point P c = [X c ,Y c ,Z c ,l] projects to point p = [x,y,l] on the image 
plane by the foil owing equations: 
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y = f Y c 

In practice, we are not able to directly access the image plane coordinates. 
I nstead we have access to an array of pixels in computer memory I n order to 
understand the relationship between image plane coordinates and the array 
of pixels, we must examine the imaging process. The image plane of a CCD 
camera is a rectangular array of discrete light sensitive elements. The output 
from each of these elements is an analog signal proportional to the amount 
of light incident upon it. The values for each of these elements are read out 
one element at a time row after row until the entire sensor array has been 
read. The analog signal is converted to digital values by internal circuitry. 
Digital cameras use this signal directly. Analog cameras require the use of 
an additional external A to D converter commonly known as a frame grabber. 
In addition, color cameras require the signals from adjacent red, green, and 
blue sensor sites to be combined. The result is an array of digital values which 
can be read into the memory of a computer. These processes introduce noise 
and ambiguity in thethe image data (eg. pixel jitter which extreme cases can 
cause the x and y axes to appear n on -orthogonal or skewed [Lenz and Tsai, 
1988]). Most camera models omit skew angle 6> xy (the angle between the x 
and y axes minus 90°). We include 9 xy but assume it is 0. A single element of 
the array in memory is commonly called a pixel. We will refer to the row and 
column number of a given pixel as y 1 and x' respectively. Several parameters 
are defined to quantify the relationship between the array in memory and the 
coordinate system of the image plane. x andy are the pixel coordinates of the 
principal point. s x and s y arethe number of pixels in memory per unit distance 
in there andy direction of the image plane. These parameters along with / and 
6» xy are intrinsic or internal camera calibration parameters. The projection of 
poi nt P c to poi nt p' = [x',y',l] in memory is described by the foil owing equations 
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or more compactly 
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C is a 3 x 3 lower triangular matrix which contains the internal camera param- 
eters. 
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The right half of FigureA-1 shows an arbitrary Cartesian coordinate system 
which we will refer to as the world or absolute coordinate system. A point 
P w in the world coordinate system is transformed into the camera centered 
coordinate system by the foil owing equation: 

p c = p w n + T* 

or 



P = P 



n o 

t* 1 



(A.2) 



Where TZ is a 3 x 3 orthonormal rotation matrix and f isalx3 translation 
vector. This matrix is commonly referred to as the extrinsic or external camera 
parameters. 

The pin-hole camera is a linear model and is not able to model nonlinear 
effects such as lens distortion. The major categories of lens distortion are: 

1. Radial distortion - the path of a light ray traveling from the object to the 
image plane through the lens is not always a straight line. 

2. Decentering distortion -the optical axis of individual lens components are 
not always col linear. 

3. Thin prism distortion -the optical axis of the lens assembly is not always 
perpendicular to the image plane. 
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Appendix B 
Gradient Derivations 



This appendix presents the symbolic gradient expressions used along with con- 
jugate gradient methods to find the optimum shift (Section B.l) and surfel nor- 
mal (Section B.2). 

B.l Gradient Expression for Updating Shifts 

For simplicity, we derive the results for a single color channel (red). To obtain 
expressions for the remaining col or channels simply substitute the appropriate 
color channel values. Equation 4.3 can be rewritten using Equations 2.9 and 
3.11 as 

e(u, v) = e T (u, v) + e g (u, v) + e h (u, v) (B.l) 



where 



*(„,„)=£ (("W, + *)-■»)' / £ !, (B.2) 



ri = r 



(&}) , (B.3) 



and 



r 2 = rftft). (B.4) 



s J 



Using linear regression to calculate n\ and d r yields 1 



x,y x,y x,y x,y ,,-, ,-% 

mr = ^^ 2 v^ v^ (B.5) 



x,y x,y x,y x,y 

and 



, x,y x,y x,y x,y /r -, ,-x 

rf r = f^ 2 * * • (B.6) 

x,2/ x,y x,y x,y 



^To simplify the presentation x,y is used in place of %Sj e S r 
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The gradient of Equation B.l is 



V UtV e(u,v) = V u , v e T (u,v) + V UjV e g (u,v) + V U)V e h (u,v) 



(B.7) 



de r de T 
du ' dv 


+ 


de g de g 
du ' dv 


+ 


de h de h 
du ' dv 



Again for simplicity, we derive only f^. To obtain an expression for ^ simply 
substitute^ for u. Differentiating Equations B.2, B.5, and B. 6 with respect tou 
gives 



&, = (Kn + jj - r 2 ) (m,fa + d -tn + f - % ) / i 



(B.8) 
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(B.9) 
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(B.10) 



u and v are shifts in the^ and yl directions respectively in imagei, therefore 
differentiating Equations B. 3 and B. 4 with respect tow produces 



(B.ll) 
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Substituting Equations B.9, B.10, B.ll, and B.12 into Equation B.8 results in 
an expression which depends only upon the red channel of the raw image data 
and the gradient of the red channel in, r 2 , and ^-). Similar results apply for 
^ as well asthegreen and blue channels. 



B.2 Gradient Expression for Updating Normals 

Again for simplicity, we derive the results for a single color channel (red). To 
obtain expressions for the remaining color channels simply substitute the ap- 
propriate color channel values. Similar to Equation B.l, Equation 4.7 can be 
rewritten as 

e(0,<f>) = e I (e,<f>) + eJ0,</>) + e h ($,<l>) 



where 



and 



E E 



<.x,y 



((m r ri + d r ) - r 2 f 
of 

r 1 = r(T(^P) + Ki 
r 2 = r(raS J; n). 






(B.13) 

(B.14) 
(B.15) 

(B.16) 




Figure B-l: Surfel with attached coordinate system. 

Figure B-l shows surfel S,- with an attached coordinate system. Points on S,- 
can be specified as an offset along xl and yl relative to P.,, 



2/Q . 



Pj + xwxl + ywyl. 



(B.17) 



w is a scale parameter which determines the spacing between points on the 
surfel. The normal is parameterized as a small rotation relative to the n,. 9 
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specifies the axis of rotation 



if = xt cos + yl sin 



(B.18) 



and 4> specify the angl e of rotati on about if. Accounti ng for the updated normal , 
Equation B.17 becomes 



l$j = Pj + xwK(9, c/))xl + yw1Z(6, (f))yl 



where 

11(9,4 



(B.19) 



(B.20) 



q x + (1 ~~ 1x) cos0 q x q y (I — cos (j)) + q z sin (j) q x q z (l — cos0) — q y sin< 



q x q y (1 - cos 0) - q z sin q£ + (1 - qj) cos 

<7x9z(l — cos 0) + q y sin g y g z (l — cos <fi) — q x sin < 

The image point to which |Sj projects is 

T(lS J ,I) = lx',y'} 
where 



q y q z (l - cos 0) + q x sin < 
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The gradient of Equation B.13 is 
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First we will examined. Differentiating Equation B.14gives 
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Differentiating Equations B.9 and B.10 and using the chain rule produces 
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The expressions for ^ and ^ are given in the last section. Differentiating 

~ ox ay -* -* 

Equations B.15 and B. 16 and using thechain ruleyields 
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and 
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Equations B.21 and B.22 produce 
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Differentiating Equation B. 19 yields 
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and 

— — = y s cos 9 - x s sin 6*. (B.33) 

Substituting Equations B.25, B.26, B.27, B.28, B.29, B.30, B.31, B.32, and B.33 
into Equation B.14 results in an expression for ^ that is easily evaluated and 
depends only upon known quantities. Specifically, the red channel of the raw 
image data (r t and r 2 ), the gradient of the red channel (f^ and ^ L ), a scale 
parameter (w), the initial surfel (P,-, a%, and yt), the best shifts (w and v), and 
the camera parameters (/, x , y , x£, y£, z£, and CO. 

The derivation for J^ is nearly identical. drR Q ( f ) istheonly significant differ- 
ence, 
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Similar results apply for the green and blue channels. 
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